---
title: "Main Analysis"
author: '26966767'
date: "2024-09-01"
output: html_document
---


```{r}
setwd("/Users/fmsaitama/Downloads/Fukumoto_APSR")

```

```{r}
library(here)
```


```{r}
library(ggplot2)
library(dplyr)
library(tidyr)
library(fixest)
library(modelsummary)
library(stargazer)
library(tidyverse)
library(igraph)
library(reshape2)
library(rgenoud)
library(glmnet)   

```


```{r}


# Reading data
# ehd: Any legislators 1936-1945: Main data used for analysis
ehd <- read.csv(here("Legislator_main.csv"))
ehd <- ehd[!(apply(ehd, 1, function(x) all(is.na(x) | x == ""))), ]

# ehe: Only those who participated in legislative sessions 1936-1942
ehe <- read.csv(here("Legislator_main.csv"))
ehe <- ehe[!(apply(ehe, 1, function(x) all(is.na(x) | x == ""))), ]
ehe <- ehe[ehe$P3642 == 1, ]
```


```{r}
# hirota: Hirota cabinet vote results
hirota <- read.csv(here("Hirota.csv"))
hirota$fr <- - hirota$Ashida
hirota <- hirota[!(apply(hirota, 1, function(x) all(is.na(x) | x == ""))), ]

# takao: Saito Takao expulsion vote result
takao <- read.csv(here("takao.csv"))
takao <- takao[!(apply(takao, 1, function(x) all(is.na(x) | x == ""))), ]

# takao and hirota: merge with ehd 
ehd <- merge(ehd, hirota, by = "name", all.x = TRUE)
ehd <- merge(ehd, takao, by = "name", all.x = TRUE)

#ehe <- merge(ehe, hirota, by = "name", all.x = TRUE)
#ehe <- merge(ehe, takao, by = "name", all.x = TRUE)

```




```{r}

# Roll call vote: Most motions are irrelevant to the paper, so each motion is not described. Specific motions are explained in the following code.

#79th diet
law79 <- read.csv(here("laws7779.csv"))
law79b <- read.csv(here("laws7779b.csv"))

#73th diet
law73 <- read.csv(here("law73.csv"))

# Sessions in 1938
law38b <- read.csv(here("law38b.csv"))

#75th Diet
law75a <- read.csv(here("law75a.csv"))
law75c <- read.csv(here("law75c.csv"))

# merge with ehd

ehd <- merge(ehd,law79,by="name",all.x=TRUE)
ehd <- merge(ehd,law79b,by="name",all.x=TRUE)
ehd <- merge(ehd,law73,by="name",all.x=TRUE)
ehd <- merge(ehd,law38b,by="name",all.x=TRUE)
ehd <- merge(ehd,law75a,by="name",all.x=TRUE)


```







```{r}
# Make a subset of legislators present in each session / vote

eh1 <- ehd[ehd$Diet36 == 1 ,]
eh1b <- ehd[ehd$Diet36 == 1 ,]
eh2 <- ehd[ehd$EDR37 == 1,]
eh2b <- ehd[ehd$B380218 ==1,]
eh2c <- ehd[ehd$B380218 ==1,]
eh3 <- ehd[ehd$Diet37 ==1,]
eh3c <- ehd[ehd$B390311 ==1,]
eh4 <- ehd[ehd$B400203 ==1,]
eh4b <- ehd[ehd$B410223 ==1,]
eh5 <- ehd[ehd$Diet37 ==1,]
eh5b <- ehd[ehd$Diet37 ==1,]
eh5c <- ehd[ehd$B411118 ==1,]
eh6 <- ehd[ehd$EDRW42 ==1,]


# Create Pro_army score of each vote explained in the manuscript for those present in such a way that they fall in 0-1 scale

eh1$Pro_army <- 1 - eh1$Rights
eh1b$Pro_army <- 1 - eh1b$Ashida
eh2$Pro_army <- eh2$Hayashi
eh2b$Pro_army <- 1 - eh2b$Prosecute38
eh2c$Pro_army <- 1 - eh2c$Anti38
eh3$Pro_army <- eh3$Faction39_Desc 
eh3c$Pro_army <- 1 - eh3c$Q39.3
eh4$Pro_army <- eh4$Takao 
eh4b$Pro_army <- 1 - eh4b$Anti41_2
eh5$Pro_army <- ifelse(eh5$Group41 == 2, 1, 0)
eh5c$Pro_army <- - eh5c$Takechi41_11_18 + 1
eh5b$Pro_army <- eh5b$Army_Grade/2
eh6$Pro_army <- ifelse(eh6$Every42 == 1, 1, 0)

# Set the date of vote 

eh1$treatment_factor <- '1937-03-22'
eh1b$treatment_factor <- '1937-03-31'
eh2$treatment_factor <- '1937-04-30'
eh2b$treatment_factor <- '1938-02-21'
eh2c$treatment_factor <- '1938-02-21'
eh3$treatment_factor <- '1939-05-30'
eh3c$treatment_factor <- '1939-03-11'
eh4$treatment_factor <- '1940-02-03'
eh4b$treatment_factor <- '1941-02-23'
eh5$treatment_factor <- '1940-10-01'
eh5c$treatment_factor <- '1941-11-18'
eh5b$treatment_factor <- '1942-01-16'
eh6$treatment_factor <- '1942-04-04'



# Create treatment period (parliamentary session) for classification

eh1$treatment_period <- -2
eh1b$treatment_period <- -1.5
eh2$treatment_period <- -1
eh2b$treatment_period <- 1
eh2c$treatment_period <- 1
eh3$treatment_period <- 1.5
eh3c$treatment_period <- 1.7
eh4$treatment_period <- 2
eh4b$treatment_period <- 3
eh5$treatment_period <- 3
eh5b$treatment_period <- 5
eh5c$treatment_period <- 4
eh6$treatment_period <- 5

# Get the long data by binding the datasets above 

# dfo is the main dataset. 

dfo <- rbind(eh1,eh2,eh2c,eh3,eh3c,eh4,eh4b,eh5,eh5c,eh6)

# dfb removes those who never sat in the parliament before 1942. Used mainly for event study graph (Incumbents)

dfb <- dfo[is.na(dfo$New42),]


# dfc: Those legislator that were present at every single event

common_names <- Reduce(intersect, list(
  eh1$name, eh2$name, eh2c$name, eh3$name, eh3c$name,
  eh4$name, eh4b$name, eh5$name, eh5c$name, eh6$name
))

# Subset dfo to only include those names
dfc<- dfo[dfo$name %in% common_names, ]



# dfg For robustness checks: Includes the extra 2 events explained in the manuscript 

dfg <- rbind(eh1,eh1b,eh2,eh2c,eh3,eh3c,eh4,eh4b,eh5,eh5c,eh6,eh5b)

# dfp: For robustness checks: Only parliamentary votes as dependent variable
dfp <- rbind(eh1,eh2c,eh3c,eh4,eh4b,eh5c)

# dff: For robustness checks: Only Factions as dependent variable
dff <- rbind(eh2,eh3,eh5,eh6) 



# Set different treatment period for sanction and procurement

dfo1 <- dfo %>% mutate(Treatment = ifelse(treatment_period >= 2, 1, 0)) # Procurement shock 1 = German Invasion of Poland
dfo2 <- dfo %>% mutate(Treatment = ifelse(treatment_period >= 3, 1, 0)) # Sanction shock 1 = First set of embargo 1940
dfo3 <- dfo %>% mutate(Treatment = ifelse(treatment_period >= 4, 1, 0)) # Sanction shock 2 = Second set of embargo + asset freeze 1941
dfo3B <- dfo3[dfo3$treatment_period!=3,]                                # The period between two sanction shocks are removed for dfo3
dfo4 <- dfo %>% mutate(Treatment = ifelse(treatment_period >= 5, 1, 0)) # Procurement shock 2 = Japanese Attack on Pearl Harbor


# Do the same for dfb, dfc, dff, dfp, dfg

dfb1 <- dfb %>% mutate(Treatment = ifelse(treatment_period >= 2, 1, 0))
dfb2 <- dfb %>% mutate(Treatment = ifelse(treatment_period >= 3, 1, 0))
dfb3 <- dfb %>% mutate(Treatment = ifelse(treatment_period >= 4, 1, 0))
dfb4 <- dfb %>% mutate(Treatment = ifelse(treatment_period >= 5, 1, 0))
dfb3B <- dfb3[dfb3$treatment_period!=3,]

dfc1 <- dfc %>% mutate(Treatment = ifelse(treatment_period >= 2, 1, 0))
dfc2 <- dfc %>% mutate(Treatment = ifelse(treatment_period >= 3, 1, 0))
dfc3 <- dfc %>% mutate(Treatment = ifelse(treatment_period >= 4, 1, 0))
dfc4 <- dfc %>% mutate(Treatment = ifelse(treatment_period >= 5, 1, 0))
dfc3B <- dfc3[dfc3$treatment_period!=3,]


dff1 <- dff %>% mutate(Treatment = ifelse(treatment_period >= 2, 1, 0))
dff2 <- dff %>% mutate(Treatment = ifelse(treatment_period >= 3, 1, 0))
dff3 <- dff %>% mutate(Treatment = ifelse(treatment_period >= 4, 1, 0))
dff4 <- dff %>% mutate(Treatment = ifelse(treatment_period >= 5, 1, 0))
dff3B <- dff3[dff3$treatment_period!=3,]

dfp1 <- dfp %>% mutate(Treatment = ifelse(treatment_period >= 2, 1, 0))
dfp2 <- dfp %>% mutate(Treatment = ifelse(treatment_period >= 3, 1, 0))
dfp3 <- dfp %>% mutate(Treatment = ifelse(treatment_period >= 4, 1, 0))
dfp4 <- dfp %>% mutate(Treatment = ifelse(treatment_period >= 5, 1, 0))
dfp3B <- dfp3[dfo3$treatment_period!=3,]

dfg1 <- dfg %>% mutate(Treatment = ifelse(treatment_period >= 2, 1, 0))
dfg2 <- dfg %>% mutate(Treatment = ifelse(treatment_period >= 3, 1, 0))
dfg3 <- dfg %>% mutate(Treatment = ifelse(treatment_period >= 4, 1, 0))
dfg4 <- dfg %>% mutate(Treatment = ifelse(treatment_period >= 5, 1, 0))
dfg3B <- dfg3[dfo3$treatment_period!=3,]



```

```{r}

# Main tables (WP Table 7 and 9) and Main graphs (WP Figure 14 and 19) are presented first, then proceeds with figures in the order they appear in paper

```



```{r}

# Table 5

### WP Table 7, Columns 1 and 2 (Sanction D-i-D)

#install.packages("fixest")

library(fixest)
# The results were produced using fixest 0.12.x, in which singletons are retained in the estimation sample.
# To reproduce the results with fixest ≥ 0.13, set fixef.rm = "perfect" (equivalent to the former infinite_coef) in all fixest models.
# This adjustment only affects the number of observations reported (222 dropped out of 4648 observations in the main specification) and does not change the estimated coefficients or standard errors.
library(modelsummary)

#  Used for Fixed Effects (and clustering SE) in every model below
head(dfo2$name) # Legislator FE
head(dfo2$treatment_factor) # Event FE 


# Run the DiD regression: Sanctioned + 1st shock: dfo2
did_model_s1 <- feols(Pro_army ~ Treatment*Sanctioned  | treatment_factor + name, data = dfo2, cluster = ~name , fixef.rm = "infinite_coef")

summary(did_model_s1)

# clustering standard errors with both individual and vote
did_model_s1b <- feols(Pro_army ~ Treatment*Sanctioned  | treatment_factor + name, data = dfo2, cluster = ~name+treatment_factor, fixef.rm = "infinite_coef")
summary(did_model_s1b)


# Run the DiD regression Sanctioned + 2nd shock: dfo3B (the period between 1st and 2nd shock removed)
did_model_s2 <- feols(Pro_army ~ Treatment*Sanctioned  | treatment_factor + name, data = dfo3B, cluster = ~name, fixef.rm = "infinite_coef")

# Summary of the model
summary(did_model_s2)

# clustering standard errors with both individual and vote
did_model_s2b <- feols(Pro_army ~ Treatment*Sanctioned  | treatment_factor + name, data = dfo3B, cluster = ~name+treatment_factor, fixef.rm = "infinite_coef")

summary(did_model_s2b)

model_list <- list("Model 1" = did_model_s1, "Model 2" = did_model_s1b, "Model 3" = did_model_s2, "Model 3" = did_model_s2b)

# Output a LaTeX table
modelsummary(model_list, output = "latex")


```


```{r}

# Table 5

# Specifications excluding those new legislators for 1942-1946 parliament (dfb): Identical results to Table 7 Columns 1 & 2

# Run the DiD regression: Sanctioned + 1st shock: dfo2
did_model_s3 <- feols(Pro_army ~ Treatment*Sanctioned  | treatment_factor + name, data = dfb2, cluster = ~name, fixef.rm = "infinite_coef")

summary(did_model_s3)

# clustering standard errors with both individual and vote
did_model_s3b <- feols(Pro_army ~ Treatment*Sanctioned  | treatment_factor + name, data = dfb2, cluster = ~name+treatment_factor, fixef.rm = "infinite_coef")
summary(did_model_s3b)


# Run the DiD regression Sanctioned + 2nd shock: dfo3B (the period between 1st and 2nd shock removed)
did_model_s4 <- feols(Pro_army ~ Treatment*Sanctioned  | treatment_factor + name, data = dfb3B, cluster = ~name, fixef.rm = "infinite_coef")

# Summary of the model
summary(did_model_s4)

# clustering standard errors with both individual and vote
did_model_s4b <- feols(Pro_army ~ Treatment*Sanctioned  | treatment_factor + name, data = dfb3B, cluster = ~name+treatment_factor, fixef.rm = "infinite_coef")

summary(did_model_s4b)

model_list <- list("Model 1" = did_model_s3, "Model 2" = did_model_s3b, "Model 3" = did_model_s4, "Model 4" = did_model_s4b)

# Output a LaTeX table
modelsummary(model_list, output = "latex")

```

```{r}

# Table 5

### WP Table 7, Columns 3 and 4 (Sanction D-i-D)

# Specifications excluding those new legislators for 1942-1946 parliament (dfb): Identical results to Table 7 Columns 1 & 2

# Run the DiD regression: Sanctioned + 1st shock: dfo2
did_model_s5 <- feols(Pro_army ~ Treatment*Sanctioned  | treatment_factor + name, data = dfc2, cluster = ~name, fixef.rm = "infinite_coef")

summary(did_model_s5)

# clustering standard errors with both individual and vote
did_model_s5b <- feols(Pro_army ~ Treatment*Sanctioned  | treatment_factor + name, data = dfc2, cluster = ~name+treatment_factor, fixef.rm = "infinite_coef")
summary(did_model_s5b)


# Run the DiD regression Sanctioned + 2nd shock: dfo3B (the period between 1st and 2nd shock removed)
did_model_s6 <- feols(Pro_army ~ Treatment*Sanctioned  | treatment_factor + name, data = dfc3B, cluster = ~name, fixef.rm = "infinite_coef")

# Summary of the model
summary(did_model_s6)

# clustering standard errors with both individual and vote
did_model_s6b <- feols(Pro_army ~ Treatment*Sanctioned  | treatment_factor + name, data = dfc3B, cluster = ~name+treatment_factor, fixef.rm = "infinite_coef")

summary(did_model_s6b)

model_list <- list("Model 1" = did_model_s5, "Model 2" = did_model_s5b, "Model 3" = did_model_s6, "Model 4" = did_model_s6b)

# Output a LaTeX table
modelsummary(model_list, output = "latex")

```

```{r}

# Table 5

### WP Table 7 in simpler form

model_list <- list("Model 1" = did_model_s1, "Model 2" = did_model_s2, "Model 3" = did_model_s5, "Model 4" = did_model_s6)

# Output a LaTeX table
modelsummary(model_list, output = "latex")

```
```{r}

### Clustered standard errors in WP Table 7 (Second row & (\num{0.055}) & (\num{0.063}) & (\num{0.053}) & (\num{0.052}) \\)

model_list <- list("Model 1" = did_model_s1b, "Model 2" = did_model_s2b, "Model 3" = did_model_s5b, "Model 4" = did_model_s6b)

# Output a LaTeX table
modelsummary(model_list, output = "latex")


```



```{r}

# Figure A1

### WP Figure 8: Sanction


did_model <- feols(Pro_army ~ Sanctioned*treatment_factor | treatment_factor + name, cluster = ~  name,  data = dfo)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)

# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Sanctioned:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE


BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)


plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
  labs(title = "Sanctioned Sectors: D-i-D Estimates with Two-way Fixed Effect",
       x = "Time Period",
       y = "Treatment Effect on Pro-Army Attitude") +
    theme_minimal()


library(stargazer)
stargazer(did_model)

print(plot)




```

```{r}



# Use dfb (incumbent only): Identical (due to fixed effects)

did_model <- feols(Pro_army ~ Sanctioned*treatment_factor | treatment_factor + name, cluster = ~  name,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)

# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Sanctioned:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE


BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)


plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
  labs(title = "Sanctioned Sectors: D-i-D Estimates with Two-way Fixed Effect",
       x = "Time Period",
       y = "Treatment Effect on Pro-Army Attitude") +
    theme_minimal()


library(stargazer)
stargazer(did_model)

print(plot)


```



```{r}

# The results that does not drop the votes between two sanctions: consistent results

# Run the DiD regression: Sanctioned + 1st shock: dfo2
did_model_s7 <- feols(Pro_army ~ Treatment*Sanctioned  | treatment_factor + name, data = dfo3, cluster = ~name)

summary(did_model_s3)

# clustering standard errors with both individual and vote
did_model_s7b <- feols(Pro_army ~ Treatment*Sanctioned  | treatment_factor + name, data = dfo3, cluster = ~name+treatment_factor)

summary(did_model_s3b)

```




```{r}

# Table 6

### WP Table 9, Columns 1 and 2 (Procurement D-i-D)

# Run the DiD regression: Procured + 1st shock (WWII)
did_model_p1 <- feols(Pro_army ~ Treatment*Procured  | treatment_factor + name, data = dfo1, cluster = ~name, fixef.rm = "infinite_coef")

summary(did_model_p1)

did_model_p1b <- feols(Pro_army ~ Treatment*Procured  | treatment_factor + name, data = dfo1, cluster = ~name+treatment_factor, fixef.rm = "infinite_coef")

summary(did_model_p1b)

# Run the DiD regression: Procured + 4th shock (Pearl Harbor)
did_model_p2 <- feols(Pro_army ~ Treatment*Procured  | treatment_factor + name, data = dfo4, cluster = ~name, fixef.rm = "infinite_coef")

summary(did_model_p2)

did_model_p2b <- feols(Pro_army ~ Treatment*Procured  | treatment_factor + name, data = dfo4, cluster = ~name+treatment_factor, fixef.rm = "infinite_coef")

summary(did_model_p2b)


model_list <- list("Model 1" = did_model_p1, "Model 2" = did_model_p1b, "Model 3" = did_model_p2, "Model 4" = did_model_p2b)

# Output a LaTeX table
modelsummary(model_list, output = "latex")

```
```{r}

# Table 6

### WP Table 9, Columns 3 and 4 (Procurement D-i-D)

# Run the DiD regression: Procured + 1st shock (WWII)
did_model_p3 <- feols(Pro_army ~ Treatment*Procured  | treatment_factor + name, data = dfc1, cluster = ~name, fixef.rm = "infinite_coef")

summary(did_model_p3)

did_model_p3b <- feols(Pro_army ~ Treatment*Procured  | treatment_factor + name, data = dfc1, cluster = ~name+treatment_factor, fixef.rm = "infinite_coef")

summary(did_model_p3b)

# Run the DiD regression: Procured + 4th shock (Pearl Harbor)
did_model_p4 <- feols(Pro_army ~ Treatment*Procured  | treatment_factor + name, data = dfc4, cluster = ~name, fixef.rm = "infinite_coef")

summary(did_model_p4)

did_model_p4b <- feols(Pro_army ~ Treatment*Procured  | treatment_factor + name, data = dfc4, cluster = ~name+treatment_factor, fixef.rm = "infinite_coef")

summary(did_model_p4b)


model_list <- list("Model 1" = did_model_p1, "Model 2" = did_model_p1b, "Model 3" = did_model_p2, "Model 4" = did_model_p2b)

# Output a LaTeX table
modelsummary(model_list, output = "latex")

```

```{r}

### WP Table 9 in simpler form

model_list <- list("Model 1" = did_model_p1, "Model 2" = did_model_p2, "Model 3" = did_model_p3, "Model 4" = did_model_p4)

# Output a LaTeX table
modelsummary(model_list, output = "latex")

```

```{r}

# Figure B1

### WP Figure 13 Procurement

BL <- c(as.Date('1937-03-22'),0,0,0,0)


did_model <- feols(Pro_army ~ Procured*treatment_factor | treatment_factor + name,cluster = ~ name,  data = dfo)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)

# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Procured:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

#BL <- c(1936.03,0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-12-08'), linetype = "dashed", color = "black") +    # Zero line for reference
          geom_vline(xintercept = as.Date('1939-09-01'), linetype = "dashed", color = "black") +    # Zero line for reference
        geom_vline(xintercept = as.Date('1937-07-07'), linetype = "dashed", color = "gray") +    # Zero line for reference
  labs(title = "Firms in 1942 Procurement list: Estimates with Two-way FE",
       x = "Time Period",
       y = "Treatment Effect (Estimate)") +
  theme_minimal()


print(plot)


```


```{r}

# WP Tables B1 and B2 (Robustness checks) follows the Figures section below

```




```{r}

### Figures



### WP Figures 1, 2: In summary section below the D-i-D graph section

### WP Figures 3 - 27: Presented in the following section

### WP Figures 7, 25: In stock market section (The last part of this markdown file)

```

```{r}

set.seed(20240301)

# Figure L1

### WP Figure 3: Network (career)

library(tidyverse)
library(igraph)
library(reshape2)

data <- read_csv(here("Legislator_main.csv"))


# Get the dummy for each prefecture of origin

data_with_dummies <- data %>%
  mutate(id = row_number()) %>% # Add an ID to preserve rows
  pivot_wider(names_from = prefecture, values_from = prefecture, values_fn = length, values_fill = 0)

data <- data_with_dummies


# Get the city council and Prefectural council dummy for each city / prefecture with more than  2 reps

data_with_dummies <- data %>%
  mutate(row_id = row_number()) %>%                 # Add a unique row ID
  pivot_wider(
    names_from = PREF_council,
    values_from = PREF_council,
    values_fn = length,
    values_fill = 0,
    names_prefix = "PC_"
  ) %>%
  pivot_wider(
    names_from = CITY_council,
    values_from = CITY_council,
    values_fn = length,
    values_fill = 0,
    names_prefix = "CC_"
  )

data <- data_with_dummies




selected_vars <- c(
  # Social / Family vairables
  "samurai"	,	"adapt"	,	"inherit"	,	

  # Universitiesand school 
  #  "univ", 
  "todai", "kyodai", "keio", "waseda", "chuo", "meiji", "hosei",
  "nichidai", "kwansei", "cadet", "hitotsubashi", "Ritsumei_Misc", "Doshisha_Misc",
  "Toyo_Misc", "Senshu_Misc", "Univ_Tohoku_Misc", "Univ_Hokkaido_Misc", "tedai",
  "shidai", "national", "science", "medical", "literature", "phd", "grad",
  "us_degree", "us_study", "uk_study", "france_study_Misc", "germany_study",
  "asia_study_Misc", "study_other", "study_abroad",  "teaching_school",
  "agri_school", "commerce_school", "crafts_school", "soldier_school",
  "other_high", "lycee",
  
  # Prefectures
  "三重", "京都", "佐賀", "兵庫", "北海道", "千葉", "和歌山", "埼玉", "大分", "大阪",
  "奈良", "宮城", "宮崎", "富山", "山口", "山形", "山梨", "岐阜", "岡山", "岩手",
  "島根", "広島", "徳島", "愛媛", "愛知", "新潟", "東京", "栃木", "沖縄", "滋賀",
  "熊本", "石川", "神奈川", "福井", "福岡", "福島", "秋田", "群馬", "茨城", "長崎",
  "長野", "青森", "静岡", "香川", "高知", "鳥取", "鹿児島",
  
  # Prefectural council
  "PC_三重県", "PC_京都府", "PC_佐賀県", "PC_兵庫県", "PC_北海道", "PC_千葉県",
  "PC_和歌山県", "PC_埼玉県", "PC_大分県", "PC_大阪府", "PC_奈良県", "PC_宮城県",
  "PC_宮崎県", "PC_富山県", "PC_山口県", "PC_山形県", "PC_山梨県", "PC_岐阜県",
  "PC_岡山県", "PC_岩手県", "PC_島根県", "PC_広島県", "PC_徳島県", "PC_愛媛県",
  "PC_愛知県", "PC_新潟県", "PC_東京府", "PC_栃木県", "PC_沖縄県", "PC_滋賀県",
  "PC_熊本県", "PC_石川県", "PC_神奈川県", "PC_福井県", "PC_福岡県", "PC_福島県",
  "PC_秋田県", "PC_群馬県", "PC_茨城県", "PC_長崎県", "PC_長野県", "PC_青森県",
  "PC_静岡県", "PC_香川県", "PC_高知県", "PC_鳥取県", "PC_鹿児島県",
    
  # City council
  "CC_下関市", "CC_久留米市", "CC_京都市", "CC_今治市", "CC_仙台市", "CC_佐世保市",
  "CC_倉敷市", "CC_八代市", "CC_八幡市", "CC_八幡浜市", "CC_函館市", "CC_別府市",
  "CC_半田市", "CC_名古屋市", "CC_呉市", "CC_和歌山市", "CC_四日市市", "CC_堺市",
  "CC_塩竈市", "CC_大分市", "CC_大垣市", "CC_大津市", "CC_大牟田市", "CC_大阪市",
  "CC_奈良市", "CC_姫路市", "CC_宇和島市", "CC_宇治山田市", "CC_宇部市",
  "CC_宇都宮市", "CC_宮古市", "CC_宮崎市", "CC_宮津市", "CC_富山市", "CC_小樽市",
  "CC_小田原市", "CC_尾道市", "CC_山口市", "CC_岡谷市", "CC_岩国市", "CC_岸和田市",
  "CC_川内市", "CC_川口市", "CC_川崎市", "CC_布施市", "CC_帯広市", "CC_平市",
  "CC_広島市", "CC_延岡市", "CC_新宮市", "CC_新潟市", "CC_日立市", "CC_旭川市",
  "CC_明石市", "CC_札幌市", "CC_東京市", "CC_松山市", "CC_松本市", "CC_松阪市",
  "CC_横浜市", "CC_横須賀市", "CC_水戸市", "CC_沼津市", "CC_津山市", "CC_津市",
  "CC_浜松市", "CC_浦和市", "CC_熊本市", "CC_神戸市", "CC_福井市", "CC_福岡市",
  "CC_福島市", "CC_米沢市", "CC_若松市", "CC_諫早市", "CC_豊橋市", "CC_足利市",
  "CC_那覇市", "CC_都城市", "CC_金沢市", "CC_鎌倉市", "CC_長岡市", "CC_長崎市",
  "CC_長野市", "CC_静岡市", "CC_首里市", "CC_高松市", "CC_高田市", "CC_高知市",
  "CC_鳥取市", "CC_鶴岡市", "CC_鹿児島市", "CC_鹿屋市",
  
   # Military variable
    "Military",   "Navy",  "General",  "Colonel",   "Lieutenant", 
  "Conscript",   "Accounting_Military",   "Clerical_military", 
  "Military_Army_Medic_Misc",  "Military_Army_Journalist_Misc",   "Military_Army_Logistics_Misc", "Military_Army_Engeneering_Misc",  "Military_Navy_Torpedo_Misc",  "Military_Army_Cannon",  "Military_Airforce_Misc",
  # "Other_Military",
  
  # Administratice office 
  "Localgov_Tokyo", "Localgov_Osaka", "Localgov_Kyoto_Misc", "Localgov_Hokkaido", 
  "Interior",  "Governor", "Interior_Central", 
  "Police", "Localgov", 
  "Secretary", "Treasury", "Taxcollect", "Monopoly_Misc", 
  "Diplomat",  "Trainministry", 
  "Teleministry", 
  "Eduministry", 
   "Statministry", 
  "Indusministry", 
  "Otherministry_Audit_Misc", "Otherministry_Army_Misc", "Otherministry_Navy_Misc", 
  "Otherministry_Reconst_Misc", "Otherministry_Legal_Misc", "Otherministry_Colonial_Misc", 
  "Otherministry_Intel_Misc", "Ministry", "Ministry_Elite", "Ministry_Counsel",
  
  # School they taught in / administered (career)
  "Primary_Principal", "Primary_Teacher", "Primary_Staff", "Primary", "Educator", 
  "Middleschool_Teacher", "Specialschool_Teacher", "Highschool_Teacher", 
  "Colonialschool_Teacher_Misc", "Girlsschool_Teacher", "Agrischool_Teacher_Misc", 
  "Commerceschool_Teacher", "Secondary", 
  "Professor", "Lecturer", "Tertiary", 
  "Tertiary_Waseda", "Tertiary_Hitotsubashi_Misc", "Tertiary_Tokyo_Misc", 
  "Tertiary_Tsukuba_Misc", "Tertiary_Keio_Misc", "Tertiary_Nichidai", 
  "Tertiary_Meiji_Misc", "Tertiary_Rikkyo_Misc", "Tertiary_Chuo_Misc", 
  "Tertiary_Takushoku_Misc", "Tertiary_Kwansei_Misc", "Tertiary_Kangaku_Misc", 
  "Tertiary_Housei", "Tertiary_Senshu_Misc", "Tertiary_Toyo_Misc", 
  "Tertiary_Daito_Misc", "Tertiary_Taisho_Misc", "Tertiary_Doshisha_Misc",
  
   # Journalist: type of job and newspapers they worked for
    "Reporter", "Chiefreporter", "Editor", "Other_Journalist", "Journalist", 
  "Yomiuri", "Tokyo_Asahi_Misc", "Osaka_Asahi_Misc", "Tokyo_Nichinichi_Misc", 
  "Osaka_Mainichi_Misc", "Asahi", "Mainichi", "Jiji", "Hochi", "Yorozu_Misc", 
  "Yamato_Misc", "Niroku_Misc", "Chugai_Misc", "Dentsu_Misc", "Chuo_Shinbun", 
  "Miyako_Shinbun_Misc", "Kokumin_Shinbun_Misc", "Localnews", 
  "Tokyo_News", "Tabloid",
 
  # Career / Occupation variables
  "Lawyer", "Prosecutor", "Judge", 
  "Doctor", 
  "Blue_Collar", "White_Collar", "Zaibatsu_Worker", "Engeneer", 
  "Noneconomic_Job", "School_Worker", "Kindergarten_Misc", 
  "Public_Sector_Job", "Shugiin_Misc", "Nichigin_Misc", "NHK_Misc", 
  "Public_Sector_Local", "Public_Sector_National", "Legal_Clark", 
  "Religious_Job", "Religious_Job_Christian_Misc", 
  "Religious_Job_Budhist_Misc", "Religious_Jon_Shintoist_Misc", 
  "Welfare_Institution", "Accountant_Misc", 
  "Writer", "Thinktank", 
  "Pharmacist_Misc", "Dentist_Misc", "Veterinarian_Misc", 
  "Manchu_Admin", "Manchu_Career", "Manchu_Railway", 
  "Korea_Admin", "Korea_Career_Misc", 
  "Taiwan_Admin", "Taiwan_Career_Misc"
  
)



# View the resulting dataset (Present in 1942)
head(data)
data <- data[!is.na(data$name)&data$EDRW42==1,]

dummy_vars <- data %>% dplyr:: select(all_of(selected_vars))

summary(dummy_vars)

# Checking if the two rows (names) share at least one '1' in any of the variables
adjacency_matrix <- (as.matrix(dummy_vars) %*% t(as.matrix(dummy_vars))) > 0

# Convert the adjacency matrix into a graph
g <- graph_from_adjacency_matrix(adjacency_matrix, mode = "undirected", diag = FALSE)

# Add names to the vertices
V(g)$name <- data$name
V(g)$party <- data$Anti42

# Define two colors for the binary party variable (0 = opposition, 1 = ruling)
party_colors <- c( "deepskyblue", "navy")  # Blue for opposition, Red for ruling (IRAA)
party_shape <- c("circle", "square") 
# Map the party variable (0, 1) to the corresponding colors
V(g)$color <- party_colors[V(g)$party + 1]  # Add +1 to match index (R uses 1-based indexing)
V(g)$shape <- party_shape[V(g)$party + 1]  # Add +1 to match index (R uses 
# summary()

adjustshape <- function(shape_vector) {
  # Map input shapes to valid igraph shapes
  sapply(shape_vector, function(shape) {
    switch(
      shape,
      "circle" = "circle",        # Use "circle" for circular nodes
      "square" = "square",        # Use "square" for square nodes
      "triangle" = "crectangle",  # Map "triangle" to filled rectangle
      "star" = "csquare",         # Map "star" to filled square
      "none" = "none",            # Invisible nodes
      "circle" # Default fallback if shape not recognized
    )
  })
}


g <- delete.vertices(g, degree(g) == 0)

layout <- layout_with_fr(g)

E(g)$color <- adjustcolor("gainsboro", , alpha.f = 0.4)

shape_legend <- c("circle", "rectangle")
color_legend <- c( "deepskyblue", "navy")
labels <- c( "Pro-Army 1942", "Anti-Army 1942")
pch_legend <- ifelse(shape_legend == "circle", 21, 22) 


# Plot the network
# plot_graph <- function() {
  plot(g, 
     layout = layout,               # Use a different layout
     vertex.size = 2,                # Size of the vertices (dots)
     vertex.shape = adjustshape(V(g)$shape),     # Shape of the vertices
     vertex.label = NA,              # Hide the labels (names) from the plot for clarity
     edge.width = 1,                 # Thickness of the edges (lines)
     edge.color = E(g)$color,            # Color of the edges
#     vertex.color = Sanctioned,          # Color of the vertices
     vertex.color = NA,        # Label color for visibility
  vertex.frame.color = adjustcolor(V(g)$color, alpha.f = 0.9), 
#     vertex.color = adjustcolor(V(g)$color, alpha.f = 0.9),
      rescale = TRUE,               # Do not automatically rescale the layout
 #      xlim = range(layout[, 1]) * 1.2,   # Adjust x-axis limits to fill space
#     ylim = range(layout[, 2]) * 1.2,   # Adjust y-axis limits to fill space
     asp = 0)                       # Aspect ratio to stretch the plot

legend(
  "topleft",                              # Position of the legend
  legend = labels,                         # Labels
  pch = pch_legend,                        # Point shapes
  pt.bg = NA,  # Background colors of the points
  col = adjustcolor(color_legend, 0.6),                           # Border color of the points
  pt.cex = 1.5,                            # Size of the legend points
  cex = 0.8,                               # Text size
  bty = "n"                                # No box around the legend
)


# Basic network metrics
cat("Number of vertices (names):", vcount(g), "\n")
cat("Number of edges (connections):", ecount(g), "\n")
cat("Is the graph connected?", is_connected(g), "\n")


```



```{r}

# Figure L2

### WP Figure 4

set.seed(20240301)
#set.seed(20250907)

library(tidyverse)
library(igraph)
library(reshape2)

data <- read_csv("Legislator_main.csv")

# 	Agridevelop_Institute_Misc	Agrifinance_Misc	Tabacco_Misc	Chiken_Misc	Fruits_Misc	Charcoal_Misc	Grain_Misc	Tea_Misc	Agriland	Cultivate	Noukai	Fisher	Forestry	Sericulture	Livestock	Agri_Interest_Any	

#Interest_Urban	Interest_Other_Misc	Interest_Wara_Misc	Interest_Pharma_Misc	Realestate_Uni_Desc	Financeuni_Desc	ElecUni_Desc	Saltuni_Desc	Fooduni_Desc	Interest_Realestate_Misc	Interest_Finance	Interest_Elec_Misc	Interest_Salt_Misc	Interest_Food_Misc	InfraUni_desc	TourismUni_Desc	MineUni_Desc	WoodUni_Desc	ConstUni_Desc	Interest_Plank	Interest_Const	Interest_Mining_Misc	Interest_Tourism_Misc	Interest_Infra_Misc	PublisherUni_Desc	EmployerUni_Desc	FoodUni_desc	ManufUni_Desc	RetailUni_Desc	TransportUni_Desc	Interest_Publisher_Misc	Interest_Employer	Interest_Food	Interest_Manuf	Interest_Commerce	Interest_Transport	TextileUni_Desc	Cotton_Uni_Misc	Silk_Uni_Misc	Linen_Uni_Misc	Wool_Uni_Misc	Layon_Uni_Misc	Leather_Uni_Misc	Silk_Any	LiquorUni_Desc	Chukin_Desc	Sankumi_desc	LiquorUni	Chukin_Misc	Sankumi	COOP_desc	CreditUni_desc	COOP	CreditUni	COC	Interest_Industry_Any	Interest_Infra_Any	Road_League_Misc	Train_Leagiue	Const_League	Interest_Social	CasteUni_Desc	StatUni_Desc	AdminUni_Desc	SchoolUni_desc	ReligionUni_Desc	WelfareUni_Desc	Union_Caste_Misc	Association_Stat_Misc	Association_Selfgovern_Misc	Interest_School_Misc	Interest_Religion_Misc	Interest_Welfare	BilateralUni_Desc	OccuUni_Desc	Interest_Bilateral_Misc	Interest_Medical_Misc	Interest_Accountant_Misc	Interest_Mayor	SportsUni_Desc	YouthUni_Desc	MayorUni_desc	PeasantUni_Desc	LaborUni_desc	Union_Sports	Union_Youth	Union_Peasants	Union_Labor	VeteranUni_Desc	DoctprUni_desc	LawyerUni_desc	Interest_Admin_Desc	Interest_Veteran	Interest_Doctor	Interest_Lawyer	Interest_Admin	Interest_Social	Price_Audit_Desc	Planning_Audit_Desc	Protection_Audit_Desc	Probation_Officer	Comissioner_Any	Commissioner_Price	Commissioner_Planning_Misc	Commissioner_Conscript_Misc	Conscript_Audit_Desc	Tax_Audit_Desc	Land_Audit_Desc	Income_Audit_Desc	Audit_Any	Audit_Tax_Misc	Audit_Land	Audit_Income	Other_Med_Desc	Commerce_Med_Desc	Debt_Med_Desc	Rental_Med_Desc	HR_Med_Desc	Tenant_Med_Desc	Mediator_Any	Mediator_Other_Misc	Mediator_Commerce_Misc	Mediator_Debt	Mediator_Rental	Mediator_HR	Mediator_Peasant	Mediator_Any	Police_Com_Desc	Fire_Com_Desc	Health_Com_Desc	Edu_Com_Desc	Community_Any	Community_Education	Community_Health_Misc	Community_Fire	Kouan_Misc	Interest_Ideological	RightUni_Desc	LeftUni_Desc	YokusanUni_Desc	Genyousha_Misc	Kouadoumei_Misc	Kokuhonsha_Misc	Daitoukyoukai_Misc	Right_Uni_Any_Misc	Interest_Yokusan_Misc	Interest_Colonial_Desc	Interest_Colonial	Interest_Colonial_Manchuria	Interest_Colonial_Emigration_Misc	Interest_Colonial_Korea_Misc	Interest_Colonial_Taiwan_Pacific_Misc		name2	job	Listed_Desc	Remainder	Other_Business	Oil_Business_Desc	Postoffice_desc	Ration_Desc	Oil_Sales	Oil_Refine	Oil_Any	Postoffice	Ration	Business_Other	Ice_Business_Desc	Pottery_Business	Colonial_Business	Family_Business	Business_Ice	Business_Pottery	Business_Colonial	Business_Family	Business_Any	Industry_Business_Desc	name2B	Procured_full	Procured	Procured_limit	Procured_Desc	Sanctioned		Industry_Other_Desc	Industry_Machine_Desc	Industry_Auto_Desc	Industry_Metal_Desc	Industry_Heavy_Desc		Industry_Other_Misc	Industry_Machine	Industry_Auto	Industry_Metal	Industry_Steel	Industry_Heavy	Industry_Precision_Misc	Industry_Electronics	Industry_Heavy_Machine_Misc	Industry_Ship_Misc	Industry_Airplane_Misc	Industry_Chemical_Desc	Industry_Paint_Desc	Fertilizer_Desc	Fertilizer	Industry_Chemical_Paint_Fertilizer_Uni	Industry_Chemical	Industry_Paint	Industry_Light_Desc	Industry_Gum_Desc	Industry_Light_Pipe_Misc	Industry_Light_Typewriter_Misc	Industry_Light_Can_Misc	Industry_Light_Bicycle_Misc	Industry_Light_Watch_Misc	Industry_Light_Other_Misc		Industry_Light	Rubber	Industry_Any	Train_Business_Desc	Shipping _Desc	Transport_Business_Desc	Bus_Desc_Desc	Transport_Train	Transport_Ship	Transport_Land	Transport_Land_Warehouse	Transport_Land_Truck	Transport_Bus	Transport_Any	Mining_Business_Desc	CoalMining_Desc	Oilfield_Desc	Mining_Coal	Mining_Oil_Misc	Mining_Other	Mining_Any	Cotton_Desc	Textile_Silk_Desc	Linen_Desc	Wool_Desc	Layon	Other_Fabric	Textile_Cotton	Textile_Silk	Textile_Silk_Likely	Textile_Cotton_Silk	Textile_Weaver	Textile_Spinner	Textile_Linen_Misc	Textile_Wool_Misc	Textile_Leather_Misc	Textile_Other_Misc	Textile_Leather_Wool_Linen_Uni	Textile_Wool_Linen	Textile_Any	Textile_Nonleather	Textile_Nonmilitary	Wholesale_Business_Desc	Retail_business_Desc	Maket_Business_Desc	Merchant_Desc	Commerce_Wholesale	Commerce_Wholesale_Trade	Commerce_Retail	Commerce_Retail_Departmentstore	Commerce_Market	Commerce_Merchant_Misc	Commerce_Any	Const_Business_Desc	Realestate_Desc	Cement_Desc	Brick_Desc	Realestate	Const_Cement	Const_Const	Const_Brick_Misc	Const_Any	Pharma_Business	Pharma	Tourism_Business_Desc	Hotel_Desc	Movie_Desc	Onsen_Desc	Gamble_Desc	Service_Tourism	Service_Hotel	Service_Movie_Misc	Service_Onsen_Misc	Service_Gamble	Service_Any		Elec_Business	Gas_Desc	Infra_Electricity	Infra_Electricity_Hydro	Infra_Electricity_Supply	Infra_Electricity_Lighting_Misc	Infra_Gas	Infra_Any	Finance_Business_Desc	Insurance_Desc	Stock_Desc	Exchange_Desc	Bank_Any	Bank_Small	Bank_Trust_Misc	Bank_Commerce_Misc	Bank_Agricultural	Bank_Saving	Bank_Mutual	Stock_Broker_Misc	Stock_Exchange_Misc	Stock_Any	Insurance_Any	Insurance_Accident_Misc	Insurance_Life	Finance_Any	Agri_Business_Desc	Tea_Business_Desc	Dairy_Business_Desc	Silk_business_Desc	Agritool_Desc	Tabacco_Desc	Food_Other	Food_Other_Rice_Misc	Food_Other_Fruit_Misc	Food_Other_Wine_Misc	Food_Other_Hokkaido_Misc	Food_Silk	Food_Farmtool_Misc	Food_Dairy	Food_Tabacco_Misc	Food_Tea_Misc	Food_Shoyu	Food_Salt_Misc	Food_Fish	Food_Starch_Misc	Food_Flour_Misc	Food_Sugar_Misc	Food_Any	Food_nonsugar	Liquor	Paper_Business_Desc	Plank_Desc	Wood_Paper	Wood_Paper_Pulp_Misc	Wood_Plank	Wood_Plank_Forest	Wood_Plank_Sawing	Wood_Plank_Charcoal_Misc	Wood_Any	Shoyu_Desc	Fishing_Business_Desc	Salt_Desc	Denpun_Desc	Flour_Desc	Liquor_Desc	Sugar_Desc	Labor_Career_Desc	Zaibatsu_Desc	Engineer_Desc	Blue_Collar_Desc	Blue_Collar	White_Collar	Zaibatsu_Worker	Engeneer	Noneconomic_Career_Desc	Schoolstaff_Desc	Kindergarten_Desc	Noneconomic_Job	School_Worker	Kindergarten_Misc	Public_Sector_Other_Desc	Public_Sector_Job	Shugiin_Misc	Nichigin_Misc	NHK_Misc	Public_Sector_Local	Public_Sector_National	Legal_Clark	Religious_Job_Desc	Religious_Job	Religious_Job_Christian_Misc	Religious_Job_Budhist_Misc	Religious_Jon_Shintoist_Misc	Welfare_Desc	Welfare_Institution	Accountant_Misc	Accountant_Career	Activist_Desc	Writer_Career	Writer	Thinktank	Ingaidan_Misc	Thinktank_Career	Ingaidan_Desc	Dentist_Pharmacist_Vet	Pharmacist_Misc	Dentist_Misc	Veterinarian_Misc	Pharmacist_Desc	Dentist_Desc	Veterinarian_Desc	Agriculture_Career	Publisher_Business	Publisher_Owner_Misc	Publisher_Executive_Misc	Publisher_Board_Misc	Publisher_Audit_Misc	Printer_Owner_Misc	Publisher_Any	Telecom_Business	Telecom_Owner_Misc	Telecom_Executive_Misc	Telecom_Board_Misc	Telecom_Advisor_Misc	Telecom_Any	Newspaper_Business	Newspaper_Owner_Misc	Newspaper_Executive_Misc	Newspaper_Board_Misc	Newspaper_Audit_Misc	Newspaper_Advisor_Misc	Newspaper_Any	School_business	School_Owner_Misc	School_Board_Misc	Univ_Board_Misc	School_Postwar_Misc	School_Any	Medical_Desc	Hospital	Doctor	Manchu_Admin_Desc	Manchu_Business_Desc	Manchu_Career_Desc	Manchu_Railway_Desc	Manchu_Admin	Manchu_Business	Manchu_Career	Manchu_Railway	Manchu_Any	Korea_Admin_Desc	Korea_Business_Desc	Korea_Career_Desc	Korea_Admin	Korea_Business	Korea_Career_Misc	Korea_Any	Taiwan_Admin_Desc	Taiwan_Admin	Taiwan_Business_Desc	Taiwan_Business_Misc	Taiwan_Career_Desc	Taiwan_Career_Misc	Taiwan_Any	Copied_Newspaper_Business	Yomiuri	Tokyo_Asahi_Misc	Osaka_Asahi_Misc	Tokyo_Nichinichi_Misc	Osaka_Mainichi_Misc	Asahi	Mainichi	Jiji	Hochi	Yorozu_Misc	Yamato_Misc	Niroku_Misc	Chugai_Misc	Dentsu_Misc	Chuo_Shinbun	Miyako_Shinbun_Misc	Kokumin_Shinbun_Misc	Localnews	Tokyo_News	Tabloid	Journalist＿Career	Reporter	Chiefreporter	Editor	Other_Journalist	Journalist	Tertiary_Career	Professor	Lecturer	Tertiary	Tertiary_Waseda	Tertiary_Hitotsubashi_Misc	Tertiary_Tokyo_Misc	Tertiary_Tsukuba_Misc	Tertiary_Keio_Misc	Tertiary_Nichidai	Tertiary_Meiji_Misc	Tertiary_Rikkyo_Misc	Tertiary_Chuo_Misc	Tertiary_Takushoku_Misc	Tertiary_Kwansei_Misc	Tertiary_Kangaku_Misc	Tertiary_Housei	Tertiary_Senshu_Misc	Tertiary_Toyo_Misc	Tertiary_Daito_Misc	Tertiary_Taisho_Misc	Tertiary_Doshisha_Misc	Secondary_Career	Middleschool_Teacher	Specialschool_Teacher	Highschool_Teacher	Colonialschool_Teacher_Misc	Girlsschool_Teacher	Agrischool_Teacher_Misc	Commerceschool_Teacher	Secondary	Primary_Career	Primary_Principal	Primary_Teacher	Primary_Staff	Primary	Educator	Law_Career	Lawyer	Prosecutor	Judge	Ministry	Ministry_Elite	Ministry_Counsel	Ministries_Career	Otherministry_Audit_Misc	Otherministry_Army_Misc	Otherministry_Navy_Misc	Otherministry_Reconst_Misc	Otherministry_Legal_Misc	Otherministry_Colonial_Misc	Otherministry_Intel_Misc	Otherministry_Misc	Otherministry_Elite_Misc	Otherministry_Counsel_Misc	Otherministry_and_Stat	Indusministry_Career	Indusministry	Indusministry_Elite_Misc	Indusministry_Counsel_Misc	Eduministry_Career	Eduministry	Eduministry_Elite_Misc	Eduministry_Counsel_Misc	Statministry_Career	Statministry	Statministry_Elite_Misc	Teleministry_Career	Teleministry	Teleministry_Elite_Misc	Teleministry_Counsel_Misc	Trainministry_Career	Trainministry	Trainministry_Elite_Misc	Trainministry_Counsel_Misc	Diplomat_Career	Diplomat	Diplomat_Elite_Misc	Diplomat_Counsel_Misc	Treasury_Career	Treasury	Taxcollect	Monopoly_Misc	Treasury_Elite	Treasury_Counsel_Misc	Secretary_Career	Secretary	Interior_Career	Police_Career	Localgov_Career	Interior	Interior_Elite	Interior_Counsel_Misc	Governor	Interior_Central	Police	Localgov	Localgov_Joyaku_Misc	Localgov_Guncho_Misc	Localgov_Shoki_Misc	Localgov_Shuji_Misc	Localgov_Shigaku_Misc	Localgov_Shunyu_Misc	Localgov_Tokyo	Localgov_Osaka	Localgov_Kyoto_Misc	Localgov_Hokkaido	Crime_Other_Desc	Crime_Election_Desc	Crime_Corruption_Desc	Crime_Violence_Desc	Crime_Fraud_Desc	Crime_Left_Desc	Crime_Right_Desc	Crime_Election	Crime_Corruption	Crime_Violence_Misc	Crime_Fraud_Misc	Crime_Left	Crime_Right	Crime_Other_Misc	Crime_Any	Military_Career	Military	Navy	General	Colonel	Lieutenant	Conscript	Accounting_Military	Clerical_military	Other_Military	Military_Army_Medic_Misc	Military_Army_Journalist_Misc	Military_Army_Logistics_Misc	Military_Army_Engeneering_Misc	Military_Navy_Torpedo_Misc	Military_Army_Cannon	Military_Airforce_Misc	army_desc	ministerial_desc	Law_Minister	Health_Minister	Commerce_Minister	Agriculture_Minister	Education_Minister	Transport_Minister	Telecommunication_Minister	Colonial_Minister	Navy_Minister	Army_Minister	Interior_Minister	Foreign_Minister	Fianance_Minister	Shidehara	CabinetBefore25	Kato	Wakatsuki1	Wakatsuki2	Tanaka	Hamaguchi	Inukai	Cabinet2532	Saito	Okada	Hirota	Cabinet3236	Hiranuma	Abe	Yonai	Konoe1	Konoe2	Cabinet3742	Tojo	Koiso	Suzuki	Cabinet_Wartime	Daijin	Daijin36	Daijin_Post_Misc	Seimu	Seimu36	Seimu_Post_Misc	Sanyo	Sanyo36	Sanyo_Post_Misc	Cabinet	Cabinet36	Cabinet_Post_Misc	ministerial2_desc	committee_desc	Speaker	Chair_Budget	Chair_Audit	Chair_Discipline	Chair_Ceremony	Chair_Petition	Committee_Chair	party_desc	Daihyo	Seicho	Kanjicho	Soumu	Party_Leader	Party_Role	collabo_desc	Yokusan_lead	Yokusan_cen	Yokusan_pref	Yokusan_gen	Purge_desc	postwar_death	Purge	Dead46	postwar		farright	extra1	extra2


selected_vars <- c(
  
  # Board membership -sector variable
 "Oil_Sales", "Oil_Refine", "Oil_Any",   "Postoffice", "Ration", "Business_Ice", "Business_Pottery",   "Business_Colonial", "Procured_full",  "Sanctioned",  "Industry_Machine", "Industry_Auto",   "Industry_Metal", "Industry_Steel", "Industry_Heavy",   "Industry_Precision_Misc", "Industry_Electronics",   "Industry_Heavy_Machine_Misc", "Industry_Ship_Misc",  "Industry_Airplane_Misc", "Fertilizer",   "Industry_Chemical",  "Industry_Paint",  "Industry_Light_Pipe_Misc", "Industry_Light_Typewriter_Misc",   "Industry_Light_Can_Misc", "Industry_Light_Bicycle_Misc",   "Industry_Light_Watch_Misc",    "Rubber", "Industry_Any",   "Transport_Train", "Transport_Ship", "Transport_Land",   "Transport_Land_Warehouse", "Transport_Land_Truck",   "Transport_Bus", "Transport_Any",   "Mining_Energy", "Mining_Oil_Misc",  "Mining_Any",   "Pharma", 
 "Textile_Cotton", "Textile_Silk", "Textile_Silk_Likely",  "Textile_Weaver", "Textile_Spinner",   "Textile_Linen_Misc", "Textile_Wool_Misc", "Textile_Leather_Misc",   "Textile_Other_Misc",  "Textile_Any",   "Commerce_Wholesale", "Commerce_Wholesale_Trade",  "Commerce_Retail", "Commerce_Retail_Departmentstore",  "Commerce_Market", "Commerce_Merchant_Misc", "Commerce_Any", "Realestate", "Const_Cement", "Const_Const", "Const_Brick_Misc", "Const_Any",  "Service_Tourism", "Service_Hotel", "Service_Movie_Misc",   "Service_Onsen_Misc", "Service_Gamble", "Service_Any",   "Infra_Electricity", "Infra_Electricity_Hydro",   "Infra_Electricity_Supply", "Infra_Electricity_Lighting_Misc",   "Infra_Fuel_Gas", "Infra_Any",   "Bank_Any", "Bank_Small", "Bank_Trust_Misc",   "Bank_Commerce_Misc", "Bank_Agricultural", "Bank_Saving",   "Bank_Mutual", "Stock_Broker_Misc", "Stock_Exchange_Misc",   "Stock_Any", "Insurance_Any", "Insurance_Accident_Misc",   "Insurance_Life", "Finance_Any" ,
 
"Food_Other_Rice_Misc", "Food_Other_Fruit_Misc",   "Food_Other_Wine_Misc", "Food_Other_Hokkaido_Misc", "Food_Silk",   "Food_Farmtool_Misc", "Food_Dairy", "Food_Tabacco_Misc",   "Food_Shoyu", "Food_Salt_Misc", "Food_Fish",   "Food_Starch_Misc", "Food_Flour_Misc", "Food_Sugar_Misc",   "Food_Any",  "Liquor",  "Wood_Paper", "Wood_Paper_Pulp_Misc", "Wood_Plank",  "Wood_Plank_Forest", "Wood_Plank_Sawing", "Wood_Plank_Charcoal_Misc",  "Wood_Any",

 #  Folllowing unspecified industries are not included
 "Industry_Other_Misc",  "Business_Family", "Mining_Other", "Food_Other", "Food_nonsugar",

  
  "Newspaper_Any",  "Telecom_Any",    "Publisher_Any",     "Hospital", "Manchu_Business", "Manchu_Railway",   "Korea_Business", "Taiwan_Business_Misc",

# The following is not included 
# "School_Any",  "School_Owner_Misc", "School_Board_Misc", "Univ_Board_Misc", 
 
  
  # Interest groups related to economic or agricultural activities
 "Otheragriculture_Misc", "Mashroom_Misc", "Agriinstitute_Other_Misc", "Agridevelop_Institute_Misc", "Agrifinance_Misc", "Tabacco_Misc",  "Chiken_Misc", "Fruits_Misc", "Charcoal_Misc", "Grain_Misc", "Tea_Misc",   "Agriland", "Cultivate", "Noukai", "Fisher", "Forestry", "Sericulture",   "Livestock", "Agri_Interest_Any", "Interest_Other_Misc",   "Interest_Wara_Misc", "Interest_Pharma_Misc", "Interest_Realestate_Misc",   "Interest_Finance", "Interest_Elec_Misc", "Interest_Salt_Misc",   "Interest_Food_Misc", "Interest_Plank", "Interest_Welfare", "Interest_Medical_Misc" ,"Union_Peasants", "Union_Labor",  "Interest_Const", "Interest_Mining_Misc", "Interest_Tourism_Misc", "Interest_Infra_Misc",   "Interest_Publisher_Misc", "Interest_Employer", "Interest_Food", "Interest_Manuf",  "Interest_Commerce", "Interest_Transport",   "Cotton_Uni_Misc", "Silk_Uni_Misc", "Linen_Uni_Misc", "Wool_Uni_Misc",  "Layon_Uni_Misc", "Leather_Uni_Misc", "Silk_Any", "LiquorUni",   "Chukin_Misc", "Sankumi", "COOP", "CreditUni", "COC",  "Interest_Industry_Any", "Interest_Infra_Any",   "Road_League_Misc", "Train_Leagiue", "Const_League"
  
 
 # The following interst groups are not included for being less relavant to economic activities
  # "Interest_Veteran", "Interest_Doctor", 
#  "Interest_Lawyer", "Interest_Admin",  "Interest_Social",  
 # "Interest_Colonial_Manchuria", "Interest_Colonial_Emigration_Misc",  "Interest_Colonial_Korea_Misc", "Interest_Colonial_Taiwan_Pacific_Misc", 
#  "Association_Stat_Misc", "Association_Selfgovern_Misc",    "Interest_Religion_Misc", 
 # "Interest_Bilateral_Misc", "Interest_Accountant_Misc", "Interest_Mayor",   "Union_Sports", "Union_Youth", 
 
 
)


# View the resulting dataset
head(data)
data <- data[!is.na(data$name)&data$EDRW42==1,]

dummy_vars <- data %>% dplyr:: select(all_of(selected_vars))

colSums(is.na(dummy_vars))

summary(dummy_vars)

# Checking if the two rows (names) share at least one '1' in any of the variables
adjacency_matrix <- (as.matrix(dummy_vars) %*% t(as.matrix(dummy_vars))) > 0

is_symmetric <- all(adjacency_matrix == t(adjacency_matrix))
print(is_symmetric)

# Convert the adjacency matrix to a graph object
g <- graph_from_adjacency_matrix(adjacency_matrix, mode = "undirected", diag = FALSE)

# Add names to the vertices
V(g)$name <- data$name
V(g)$party <- data$Anti42

# Define two colors for the binary party variable (0 = opposition, 1 = ruling)
party_colors <- c( "deepskyblue", "navy")  # Blue for opposition, Red for ruling
party_shape <- c("circle", "square") 
# Map the party variable (0, 1) to the corresponding colors
V(g)$color <- party_colors[V(g)$party + 1]  # Add +1 to match index (R uses 1-based indexing)
V(g)$shape <- party_shape[V(g)$party + 1]  # Add +1 to match index (R uses 
# summary()


adjustshape <- function(shape_vector) {
  # Map input shapes to valid igraph shapes
  sapply(shape_vector, function(shape) {
    switch(
      shape,
      "circle" = "circle",        # Use "circle" for circular nodes
      "square" = "square",        # Use "square" for square nodes
      "triangle" = "crectangle",  # Map "triangle" to filled rectangle
      "star" = "csquare",         # Map "star" to filled square
      "none" = "none",            # Invisible nodes
      "circle" # Default fallback if shape not recognized
    )
  })
}


g <- delete.vertices(g, degree(g) == 0)

layout <- layout_with_fr(g)

E(g)$color <- adjustcolor("gainsboro", , alpha.f = 0.4)

shape_legend <- c("circle", "rectangle")
color_legend <- c( "deepskyblue", "navy")
labels <- c( "Pro-Army 1942", "Anti-Army 1942")
pch_legend <- ifelse(shape_legend == "circle", 21, 22) 


# Plot the network
plot(g, 
     layout = layout,               # Use a different layout
     vertex.size = 2,                # Size of the vertices (dots)
     vertex.shape = adjustshape(V(g)$shape),     # Shape of the vertices
     vertex.label = NA,              # Hide the labels (names) from the plot for clarity
     edge.width = 1,                 # Thickness of the edges (lines)
     edge.color = E(g)$color,            # Color of the edges
#     vertex.color = Sanctioned,          # Color of the vertices
     vertex.color = NA,        # Label color for visibility
  vertex.frame.color = adjustcolor(V(g)$color, alpha.f = 0.9), 
#     vertex.color = adjustcolor(V(g)$color, alpha.f = 0.9),
      rescale = TRUE,               # Do not automatically rescale the layout
 #      xlim = range(layout[, 1]) * 1.2,   # Adjust x-axis limits to fill space
#     ylim = range(layout[, 2]) * 1.2,   # Adjust y-axis limits to fill space
     asp = 0)                       # Aspect ratio to stretch the plot

legend(
  "topleft",                              # Position of the legend
  legend = labels,                         # Labels
  pch = pch_legend,                        # Point shapes
  pt.bg = NA,  # Background colors of the points
  col = adjustcolor(color_legend, 0.6),                           # Border color of the points
  pt.cex = 1.5,                            # Size of the legend points
  cex = 0.8,                               # Text size
  bty = "n"                                # No box around the legend
)

# Basic network metrics
cat("Number of vertices (names):", vcount(g), "\n")
cat("Number of edges (connections):", ecount(g), "\n")
cat("Is the graph connected?", is_connected(g), "\n")


```





```{r}

# Mean pro-army (Hayashi) score in 1937

group_means <- eh2 %>%
  summarise(
    Overall = mean(Hayashi == 1),
    No_Business_Ties = mean(Hayashi[Business_Any == 0] == 1),
    Infrastructure = mean(Hayashi[Infra_Any == 1] == 1),
    Industry = mean(Hayashi[Industry_Any == 1] == 1),
    Transport = mean(Hayashi[Transport_Any == 1] == 1),
#    Textile = mean(Hayashi[Textile_Any == 1] == 1),
    Service = mean(Hayashi[Service_Any == 1] == 1),
    Food  = mean(Hayashi[Food_Any == 1] == 1),
    Newspaper = mean(Hayashi[Newspaper_Any == 1] == 1),
    Procured = mean(Hayashi[Procured == 1] == 1),
    Sanctioned = mean(Hayashi[Sanctioned == 1] == 1),
    Finance = mean(Hayashi[Finance_Any == 1] == 1),
    Business_Executives = mean(Hayashi[Business_Any == 1] == 1)
  ) %>%
  pivot_longer(cols = everything(), names_to = "Group", values_to = "Proportion")

group_means$Group <- factor(group_means$Group, levels = group_means$Group[order(group_means$Proportion)])


# Bar plot of proportions
plot<- ggplot(group_means, aes(y = Group, x = Proportion, fill = Group)) +
  geom_bar(stat = "identity", color = "black", show.legend = FALSE)  +
  geom_vline(xintercept = mean(eh2$Hayashi == 1), linetype = "dashed", color = "gray") +
  labs(
    title = "",
    x = "",
    y = ""
  ) +
  scale_fill_manual(values = c("white", "gray", "gray", "white","white", "white", "white", "lightgray", "white", "white","white","white")) +
  theme_minimal()

print(plot)


```

```{r}

# Mean pro-army (Yokusan) score in 1942 (Among those who sat in 1937-1942 Diet sessions)

eh7 <- eh6[eh6$Diet37 ==1,]

group_means <- eh7 %>%
  summarise(
    Overall = mean(Every42 == 1),
    No_Business_Ties = mean(Every42[Business_Any == 0] == 1),
    Infrastructure = mean(Every42[Infra_Any == 1] == 1),
    Industry = mean(Every42[Industry_Any == 1] == 1),
    Transport = mean(Every42[Transport_Any == 1] == 1),
#    Textile = mean(Every42[Textile_Any == 1] == 1),
    Service = mean(Every42[Service_Any == 1] == 1),
    Food  = mean(Every42[Food_Any == 1] == 1),
    Newspaper = mean(Every42[Newspaper_Any == 1] == 1),
    Procured = mean(Every42[Procured == 1] == 1),
    Sanctioned = mean(Every42[Sanctioned == 1] == 1),
    Finance = mean(Every42[Finance_Any == 1] == 1),
    Business_Executives = mean(Every42[Business_Any == 1] == 1)
  ) %>%
  pivot_longer(cols = everything(), names_to = "Group", values_to = "Proportion")

group_means$Group <- factor(group_means$Group, levels = group_means$Group[order(group_means$Proportion)])


# Bar plot of proportions
plot<- ggplot(group_means, aes(y = Group, x = Proportion, fill = Group)) +
  geom_bar(stat = "identity", color = "black", show.legend = FALSE)  +
  geom_vline(xintercept = mean(eh7$Every42 == 1), linetype = "dashed", color = "gray") +
  labs(
    title = "",
    x = "",
    y = ""
  ) +
  scale_fill_manual(values = c("gray", "white", "white", "white","lightgray", "white", "white", "white", "white", "white","white","gray")) +
  theme_minimal()

print(plot)


```




```{r}

# Figure 1

### WP Figure 5: Business-Related Variables: T-Test Results for 1937 General Hayashi Coalition
# Among those who were in the Diet 1937-1942

eh8 <- eh2[eh2$Diet37==1,]

t_test_results <- lapply(c("Sanctioned", "Procured", "Business_Any", "Infra_Any", 
                           "Industry_Any", "Transport_Any", "Service_Any", 
                           "Food_Any", "Newspaper_Any", "Finance_Any"), function(var) {
  t_test <- t.test(
    x = eh8$Hayashi[eh8[[var]] == 1],
    y = eh8$Hayashi[eh8[[var]] == 0],
    alternative = "two.sided"
  )
  cleaned_var <- gsub("_Any", "", var)
  data.frame(
    Variable = cleaned_var,
    t_statistic = t_test$statistic,
    p_value = t_test$p.value,
    mean_diff = diff(t_test$estimate)
  )
})


# Combine results into a single data frame
t_test_results_df <- do.call(rbind, t_test_results)

t_test_results_df$Variable
t_test_results_df$t_statistic


# Visualize

plot <- ggplot(t_test_results_df, aes(x = reorder(Variable, t_statistic), y = t_statistic, fill = Variable)) +
  geom_bar(stat = "identity") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
  coord_flip() +
  labs(
#    title = "T-Test Results for Business x Pro-Army Factions: 1937 Election",
    x = "",
    y = "T-Statistic"
  ) +
    scale_fill_manual(values = c("gray","lightgray","lightgray","lightgray","lightgray", "lightgray","darkgray","darkgray","lightgray","lightgray")) +
  theme_minimal() +
  theme(legend.position = "none")  # Removes the legend


print(plot)


```

```{r}

# The same t-test to WP Figure 5, but including those who lost seat in 1937 election

t_test_results <- lapply(c("Sanctioned", "Procured", "Business_Any", "Infra_Any", 
                           "Industry_Any", "Transport_Any", "Service_Any", 
                           "Food_Any", "Newspaper_Any", "Finance_Any"), function(var) {
  t_test <- t.test(
    x = eh2$Hayashi[eh2[[var]] == 1],
    y = eh2$Hayashi[eh2[[var]] == 0],
    alternative = "two.sided"
  )
  data.frame(
    Variable = var,
    t_statistic = t_test$statistic,
    p_value = t_test$p.value,
    mean_diff = diff(t_test$estimate)
  )
})


# Combine results into a single data frame
t_test_results_df <- do.call(rbind, t_test_results)

t_test_results_df$Variable
t_test_results_df$t_statistic



# Create a visualization using ggplot2

ggplot(t_test_results_df, aes(x = reorder(Variable, t_statistic), y = t_statistic, fill = Variable)) +
  geom_bar(stat = "identity") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
  coord_flip() +
  labs(
    title = "T-Test Results for Variables vs Hayashi",
    x = "",
    y = "T-Statistic"
  ) +
    scale_fill_manual(values = c("gray","lightgray","lightgray","lightgray","lightgray", "lightgray","darkgray","darkgray","lightgray","lightgray")) +
  theme_minimal() +
  theme(legend.position = "none")  # Removes the legend
```


```{r}

# Figure 2

### WP Figure 6: Business-Related Variables: T-Test Results for 1942 Army Endorsement

t_test_results <- lapply(c("Sanctioned", "Procured", "Business_Any", "Infra_Any", 
                           "Industry_Any", "Transport_Any", "Service_Any", 
                           "Food_Any", "Newspaper_Any", "Finance_Any"), function(var) {
  t_test <- t.test(
    x = eh7$Every42[eh7[[var]] == 1],
    y = eh7$Every42[eh7[[var]] == 0],
    alternative = "two.sided"
  )
  cleaned_var <- gsub("_Any", "", var)
  data.frame(
    Variable = cleaned_var,
    t_statistic = t_test$statistic,
    p_value = t_test$p.value,
    mean_diff = diff(t_test$estimate)
  )
})


# Combine results into a single data frame
t_test_results_df <- do.call(rbind, t_test_results)

t_test_results_df$Variable
t_test_results_df$t_statistic



# Create a visualization using ggplot2

plot <- ggplot(t_test_results_df, aes(x = reorder(Variable, t_statistic), y = t_statistic, fill = Variable)) +
  geom_bar(stat = "identity") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
  coord_flip() +
  labs(
#    title = "T-Test Results for Business x Pro-Army Factions: 1942 Election",
    x = "",
    y = "T-Statistic"
  ) +
    scale_fill_manual(values = c("gray","lightgray","lightgray","lightgray","lightgray", "lightgray","darkgray","darkgray","lightgray","lightgray")) +
  theme_minimal() +
  theme(legend.position = "none")  # Removes the legend


print(plot)



```



```{r}

# The same t-test to WP Figure 6, but including those who won seat for the first time in 1942 election

t_test_results <- lapply(c("Sanctioned", "Procured", "Business_Any", "Infra_Any", 
                           "Industry_Any", "Transport_Any", "Service_Any", 
                           "Food_Any", "Newspaper_Any", "Finance_Any"), function(var) {
  t_test <- t.test(
    x = eh6$Every42[eh6[[var]] == 1],
    y = eh6$Every42[eh6[[var]] == 0],
    alternative = "two.sided"
  )
  data.frame(
    Variable = var,
    t_statistic = t_test$statistic,
    p_value = t_test$p.value,
    mean_diff = diff(t_test$estimate)
  )
})


# Combine results into a single data frame
t_test_results_df <- do.call(rbind, t_test_results)

t_test_results_df$Variable
t_test_results_df$t_statistic



# Create a visualization using ggplot2

ggplot(t_test_results_df, aes(x = reorder(Variable, t_statistic), y = t_statistic)) +
  geom_bar(stat = "identity") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
  coord_flip() +
  labs(
    title = "",
    x = "Variable",
    y = "T-Statistic"
  ) +
    scale_fill_manual(values = c("gray", "white", "white", "white","lightgray", "white", "white", "white", "white", "white")) +
  theme_minimal()

```



```{r}

t_test <- t.test(
  x = eh2$Hayashi[eh2$Sanctioned == 1],
  y = eh2$Hayashi[eh2$Sanctioned == 0],
  alternative = "two.sided"
)
print(t_test)

```




```{r}

# Average Pro-Army score per event

average_data <- dfo %>%
  group_by(treatment_factor) %>%
  summarise(Average_Pro_Army = sum(Pro_army, na.rm = TRUE))

average_data

average_data <- dfo %>%
  group_by(treatment_factor) %>%
  summarise(Average_Pro_Army = sum(1-Pro_army, na.rm = TRUE))

average_data

```




```{r}

# Custom Graph Setting

custom_theme <- theme_minimal() +
  theme(
    plot.title = element_text(size = 20),  # Larger title
    axis.text = element_text(size = 14),                 # Adjust axis text size
    axis.title = element_text(size = 14),                # Adjust axis title size
    panel.grid = element_line(size = 0.5),               # Thicker grid lines
    panel.grid.major = element_line(size = 0.5, color = "grey"), # Adjust major grid line thickness
    panel.grid.minor = element_line(size = 0.25, color = "lightgrey") # Adjust minor grid line thickness
  )


```




```{r}

### Figure O1 / WP Figure 7 is in stock market section below

```






```{r}

# Figure A1

### WP Figure 8: Sanction


did_model <- feols(Pro_army ~ Sanctioned*treatment_factor | treatment_factor + name, cluster = ~  name,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)

# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Sanctioned:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE


BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)


plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
  labs(title = "Sanctioned Sectors: D-i-D Estimates with Two-way Fixed Effect",
       x = "Time Period",
       y = "Treatment Effect on Pro-Army Attitude") +
    theme_minimal()


library(stargazer)
stargazer(did_model)

print(plot)




```


```{r}

# Figure A2A

### WP Figure 9A Petrochemical


did_model <- feols(Pro_army ~ Industry_Chemical*treatment_factor | treatment_factor + name,cluster = ~   name,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)


# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Industry_Chemical:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
#       geom_vline(xintercept = 1941.5, linetype = "dashed", color = "blue") +    # Zero line for reference
  labs(title = "Petrochemical",
       x = "",
       y = "") +
    custom_theme

print(plot)


```

```{r}

# Figure A2B

# WP Figure 9B Steel

did_model <- feols(Pro_army ~ Industry_Steel*treatment_factor | treatment_factor + name, cluster = ~  name, data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)

# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Industry_Steel:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
  labs(title = "Steel",
       x = "",
       y = "") +
  custom_theme

print(plot)




```

```{r}

# Figure A2C

### WP Figure 9C Cotton Textile


did_model <- feols(Pro_army ~ Textile_Cotton*treatment_factor | treatment_factor + name,cluster = ~   name,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)

coef_names

# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Textile_Cotton:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
#       geom_vline(xintercept = 1941.5, linetype = "dashed", color = "blue") +    # Zero line for reference
 labs(title = "Cotton Textile",
       x = "",
       y = "") +
    custom_theme

print(plot)



```

```{r}

# Figure A2D

### WP Figure 9D Textile (Any)


did_model <- feols(Pro_army ~ Textile_Any*treatment_factor | treatment_factor + name,cluster = ~   name,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)


# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Textile_Any:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
#       geom_vline(xintercept = 1941.5, linetype = "dashed", color = "blue") +    # Zero line for reference
  labs(title = "Textile (Any)",
       x = "",
       y = "") +
    custom_theme

print(plot)



```
```{r}

# Figure A2E

### WP Figure 9E Stock broaker


did_model <- feols(Pro_army ~ Stock_Broker_Misc*treatment_factor | treatment_factor + name,cluster = ~   name,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)


# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Stock_Broker_Misc:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
#       geom_vline(xintercept = 1941.5, linetype = "dashed", color = "blue") +    # Zero line for reference
  labs(title = "Stock Broker",
       x = "",
       y = "") +
    custom_theme

print(plot)



```

```{r}

# Figure A2F

### WP Figure 9F International Trade


did_model <- feols(Pro_army ~ Commerce_Wholesale_Trade*treatment_factor | treatment_factor + name,cluster = ~   name,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)


# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Commerce_Wholesale_Trade:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
#       geom_vline(xintercept = 1941.5, linetype = "dashed", color = "blue") +    # Zero line for reference
  labs(title = "International Trade",
       x = "",
       y = "") +
    custom_theme

print(plot)


```

```{r}

# Figure A2G

### WP Figure 9G Petroleun (Any)


did_model <- feols(Pro_army ~ Oil_Any*treatment_factor | treatment_factor + name,cluster = ~   name,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)


# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Oil_Any:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
#       geom_vline(xintercept = 1941.5, linetype = "dashed", color = "blue") +    # Zero line for reference
  labs(title = "Petroleum (refiners and sellers)",
       x = "",
       y = "") +
   custom_theme

print(plot)



```

```{r}

# Figure A2H

### WP Figure 9H Oil seller vs Oil Refiner Comparison

did_model <- feols(Pro_army ~  Oil_Sales*Oil_Refine*treatment_factor  | treatment_factor+name,cluster = ~  name,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)


# Create a data frame for plotting
plot_data1 <- data.frame(
  Time = ifelse(grepl("^Oil_Sales:treatment_factor", coef_names),
                          as.Date(gsub("Oil_Sales:treatment_factor", "", coef_names)),
                          NA),
  Estimate = coefficients,
  SE = se
)


# Calculate confidence intervals (95%)
plot_data1$CI_lower <- plot_data1$Estimate - 1.96 * plot_data1$SE
plot_data1$CI_upper <- plot_data1$Estimate + 1.96 * plot_data1$SE
#BL <- c(1936.03,0,0,0,0)
plot_data1 <- rbind(plot_data1,BL)
plot_data1$Type <- 'Oil Sellers'

plot_data2 <- data.frame(
  Time = ifelse(grepl("^Oil_Refine:treatment_factor", coef_names),
                          as.Date(gsub("Oil_Refine:treatment_factor", "", coef_names), format = "%Y-%m-%d"),
                          as.Date(NA)),
  Estimate = coefficients,
  SE = se
)


# Calculate confidence intervals (95%)
plot_data2$CI_lower <- plot_data2$Estimate - 1.96 * plot_data2$SE
plot_data2$CI_upper <- plot_data2$Estimate + 1.96 * plot_data2$SE
#BL <- c(1936.03,0,0,0,0)
plot_data2 <- rbind(plot_data2,BL)
plot_data2$Type <- 'Oil Refiners'

plot_data <- rbind(plot_data1,plot_data2)

library(ggplot2)

plot_data <- plot_data[!is.na(plot_data$Time), ]
plot_data$Time <- as.Date(plot_data$Time)

dodge_width <- 30

plot <- ggplot(plot_data, aes(x = Time, y = Estimate, color = Type, shape = Type)) +
  geom_point(position = position_dodge(width = dodge_width) )+                    # Points for the estimates
  geom_line(aes(linetype = Type),  position = position_dodge(width = dodge_width))+                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2,size=0.8,   position = position_dodge(width = dodge_width)) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
  labs(title = "Oil Sellers vs Oil Refiner",
       x = "",
       y = "") +
  scale_color_manual(values = c( "Oil Refiners" = "grey", "Oil Sellers" = "black"), breaks = c("Oil Refiners", "Oil Sellers")) +   # Reverse order
#  scale_shape_manual(values = c("Procured" = 16, "Sanctioned" = 17), breaks = c("Sanctioned", "Procured")) +
#  scale_linetype_manual(values = c("Oil Sellers" = "solid", "Oil Refiners" = "dashed"), breaks = c("Oil Sellers", "Oil Refiners")) 
    custom_theme  +
  theme(legend.position = "bottom",
       # Enlarge legend title
    legend.text = element_text(size = 16),                  # Enlarge legend labels
    legend.key.size = unit(1.5, "lines"),                   # Increase legend key size
    legend.key.width = unit(1, "cm")                        # Adjust legend key width
  )
#scale_color_manual(values = c("blue", "green")) +  # Custom colors for cotton and linen
#scale_shape_manual(values = c(16, 17))             # Custom shapes for cotton and linen


print(plot)


```

```{r}

# Oil sellers only: Comparison with refiner reported in Figure 15H


did_model <- feols(Pro_army ~ Oil_Sales*treatment_factor | treatment_factor + name,cluster = ~   name,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)


# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Oil_Sales:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
#       geom_vline(xintercept = 1941.5, linetype = "dashed", color = "blue") +    # Zero line for reference
  labs(title = "Petroleum Sellers (Excluding Refiners)",
       x = "",
       y = "") +
    custom_theme

print(plot)



```

```{r}

# Textile spinners only: not reported in the paper, but included in Textile (Any) Figue 15B


did_model <- feols(Pro_army ~ Textile_Spinner*treatment_factor | treatment_factor + name,cluster = ~   name,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)

# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Textile_Spinner:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
#       geom_vline(xintercept = 1941.5, linetype = "dashed", color = "blue") +    # Zero line for reference
  labs(title = "Textile: Spinner (mostly silk)",
       x = "",
       y = "") +
   custom_theme

print(plot)



```

```{r}

# Figure A3A

### WP Figure 10A Insurance

did_model <- feols(Pro_army ~ Insurance_Any*treatment_factor | treatment_factor + name,cluster = ~ name ,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)

# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Insurance_Any:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference  labs(title = "Difference-in-Differences Estimates Over Time",
  labs(title = "Insurance",
       x = "",
       y = "") +
  custom_theme

print(plot)


```


```{r}

# Figure A3B

### WP Figure 10B Service sector


did_model <- feols(Pro_army ~ Service_Any*treatment_factor | treatment_factor + name,cluster = ~ name ,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)

# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Service_Any:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference  labs(title = "Difference-in-Differences Estimates Over Time",
  labs(title = "Service Sector",
       x = "",
       y = "") +
  custom_theme

print(plot)




```


```{r}

# Figure A3C

### WP Figure 10C Retail


did_model <- feols(Pro_army ~ Commerce_Retail*treatment_factor | treatment_factor + name,cluster = ~ name ,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)

# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Commerce_Retail:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference  labs(title = "Difference-in-Differences Estimates Over Time",
  labs(title = "Retail",
       x = "",
       y = "") +
  custom_theme

print(plot)


```



```{r}

# Figure A3D

### WP Figure 10D Construction


did_model <- feols(Pro_army ~ Const_Any*treatment_factor | treatment_factor + name,cluster = ~   name,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)


# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Const_Any:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
#       geom_vline(xintercept = 1941.5, linetype = "dashed", color = "blue") +    # Zero line for reference
  labs(title = "Construction",
       x = "",
       y = "") +
   custom_theme

print(plot)




```

```{r}


# Figure A3E

### WP Figure 10E: Railroad


did_model <- feols(Pro_army ~ Transport_Train*treatment_factor | treatment_factor + name,cluster = ~   name,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)


# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Transport_Train:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
#       geom_vline(xintercept = 1941.5, linetype = "dashed", color = "blue") +    # Zero line for reference
  labs(title = "Railroad",
       x = "",
       y = "") +
  custom_theme

print(plot)



```

```{r}

# Figure A3F

### WP Figure 10F: Truck and Bus


did_model <- feols(Pro_army ~ Transport_Track_and_Bus*treatment_factor | treatment_factor + name,cluster = ~   name,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)

coef_names

# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Transport_Track_and_Bus:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
#       geom_vline(xintercept = 1941.5, linetype = "dashed", color = "blue") +    # Zero line for reference
 labs(title = "Truck and Bus",
       x = "",
       y = "") +
    custom_theme

print(plot)



```

```{r}

# Figure A3G

### WP Figure 10G Newspaper

did_model <- feols(Pro_army ~ Newspaper_Any*treatment_factor | treatment_factor + name,cluster = ~   name,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)


# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Newspaper_Any:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
#       geom_vline(xintercept = 1941.5, linetype = "dashed", color = "blue") +    # Zero line for reference
  labs(title = "Newspaper",
       x = "",
       y = "") +
   custom_theme

print(plot)




```

```{r}

# Figure A3H

### WP Figure 10H: Publishing industry

did_model <- feols(Pro_army ~ Publisher_Any*treatment_factor | name  + treatment_factor,cluster = ~ name ,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)

coef_names

# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Publisher_Any:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference  labs(title = "Difference-in-Differences Estimates Over Time",
  labs(title = "Publisher",
       x = "",
       y = "") +
  custom_theme


print(plot)


```


```{r}

# Commerce: Not reported in paper (no effect, retail (subset) is repoted )


did_model <- feols(Pro_army ~ Commerce_Any*treatment_factor | treatment_factor + name,cluster = ~ name ,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)

# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Commerce_Any:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference  labs(title = "Difference-in-Differences Estimates Over Time",
  labs(title = "Commerce (Any)",
       x = "",
       y = "") +
  custom_theme

print(plot)



```

```{r}

# Wholesale industry (Oroshi): Not reported in paper (no effect, difficult to know sanction exposure)


did_model <- feols(Pro_army ~ Commerce_Wholesale*treatment_factor | treatment_factor + name,cluster = ~ name ,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)

# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Commerce_Wholesale:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference  labs(title = "Difference-in-Differences Estimates Over Time",
  labs(title = "Wholesale",
       x = "",
       y = "") +
  custom_theme

print(plot)




```


```{r}

# Real estate: not reported in paper, no effect


did_model <- feols(Pro_army ~ Realestate*treatment_factor | treatment_factor + name,cluster = ~   name,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)


# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Realestate:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
#       geom_vline(xintercept = 1941.5, linetype = "dashed", color = "blue") +    # Zero line for reference
  labs(title = "Real Estate",
       x = "",
       y = "") +
    custom_theme

print(plot)



```


```{r}

# Railway - warehouse vs Truck - Bus horse race: not reported: (neither is signicant)

BL <- c(as.Date('1937-03-22'),0,0,0,0)

did_model <- feols(Pro_army ~  Transport_Train_and_Warehouse*Transport_Nonship*treatment_factor  | treatment_factor+name,cluster = ~  name,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)


# Create a data frame for plotting
plot_data1 <- data.frame(
   Time = ifelse(grepl("^Transport_Train_and_Warehouse:treatment_factor", coef_names),
                          as.Date(gsub("Transport_Train_and_Warehouse:treatment_factor", "", coef_names), format = "%Y-%m-%d"),
                          as.Date(NA)),
#  Time = as.numeric(gsub("Textile_Wool_Linen:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)


# Calculate confidence intervals (95%)
plot_data1$CI_lower <- plot_data1$Estimate - 1.96 * plot_data1$SE
plot_data1$CI_upper <- plot_data1$Estimate + 1.96 * plot_data1$SE

plot_data1 <- rbind(plot_data1,BL)
plot_data1$Type <- 'Railroad and Warehouse'

plot_data2 <- data.frame(
    Time = ifelse(grepl("^Transport_Nonship:treatment_factor", coef_names),
                          as.Date(gsub("Transport_Nonship:treatment_factor", "", coef_names), format = "%Y-%m-%d"),
                          as.Date(NA)),
  Estimate = coefficients,
  SE = se
)

plot_data2

# Calculate confidence intervals (95%)
plot_data2$CI_lower <- plot_data2$Estimate - 1.96 * plot_data2$SE
plot_data2$CI_upper <- plot_data2$Estimate + 1.96 * plot_data2$SE

plot_data2 <- rbind(plot_data2,BL)
plot_data2$Type <- 'Truck and Bus'



plot_data <- rbind(plot_data1,plot_data2)

library(ggplot2)

plot_data <- plot_data[!is.na(plot_data$Time), ]
plot_data$Time <- as.Date(plot_data$Time)

dodge_width <- 30

plot <- ggplot(plot_data, aes(x = Time, y = Estimate, color = Type, shape = Type)) +
  geom_point(position = position_dodge(width = dodge_width)) +                    # Points for the estimates
  geom_line(aes(linetype = Type), position = position_dodge(width = dodge_width)) + # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, position = position_dodge(width = dodge_width)) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size = 2.0) +    # Zero line for reference
  geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Vertical reference line
  geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Vertical reference line
  labs(
    title = "Railroad-Warehouse vs Truck-Bus",
    x = "",
    y = ""
  ) +
  scale_color_manual(
    values = c("Railroad and Warehouse" = "black", "Truck and Bus" = "darkgrey"),
    breaks = c("Railroad and Warehouse", "Truck and Bus")
  ) +  
  scale_shape_manual(
    values = c("Railroad and Warehouse" = 16, "Truck and Bus" = 17),
    breaks = c("Railroad and Warehouse", "Truck and Bus")
  ) +  # Define shape for each type
  scale_linetype_manual(
    values = c("Railroad and Warehouse" = "solid", "Truck and Bus" = "dashed"),
    breaks = c("Railroad and Warehouse", "Truck and Bus")
  ) +  # Define line type for each type
  custom_theme  +
  theme(legend.position = "bottom",
       # Enlarge legend title
    legend.text = element_text(size = 16),                  # Enlarge legend labels
    legend.key.size = unit(1.5, "lines"),                   # Increase legend key size
    legend.key.width = unit(1, "cm")                        # Adjust legend key width
  )


print(plot)





```




```{r}

# Train and Warehouse (often combined): Not reported, but Train (Railways) is in Figure 16

did_model <- feols(Pro_army ~ Transport_Train_and_Warehouse*treatment_factor | treatment_factor + name,cluster = ~   name,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)

coef_names

# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Transport_Train_and_Warehouse:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
#       geom_vline(xintercept = 1941.5, linetype = "dashed", color = "blue") +    # Zero line for reference
 labs(title = "Train and Warehouse",
       x = "",
       y = "") +
    custom_theme

print(plot)





```





```{r}

# Land transport: not reported in the paper (no effect, no prior), but Transpoert is in Figure 16


did_model <- feols(Pro_army ~ Transport_Nonship*treatment_factor | treatment_factor + name,cluster = ~   name,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)

coef_names

# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Transport_Nonship:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
#       geom_vline(xintercept = 1941.5, linetype = "dashed", color = "blue") +    # Zero line for reference
 labs(title = "Land Transport",
       x = "",
       y = "") +
    custom_theme

print(plot)




```

```{r}

# Figure A4A

### WP Figure 11A Paper and pulp industry

did_model <- feols(Pro_army ~ Wood_Paper*treatment_factor | name  + treatment_factor,cluster = ~ name ,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)

coef_names

# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Wood_Paper:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference  labs(title = "Difference-in-Differences Estimates Over Time",
  labs(title = "Paper",
       x = "",
       y = "") +
  custom_theme


print(plot)


```

```{r}

# Figure A4B

### WP Figure 11B Forestrty industry


did_model <- feols(Pro_army ~ Fertilizer*treatment_factor | name  + treatment_factor,cluster = ~ name ,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)

coef_names

# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Fertilizer:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference  labs(title = "Difference-in-Differences Estimates Over Time",
  labs(title = "Fertilizer",
       x = "",
       y = "") +
  custom_theme


print(plot)



```

```{r}

# Figure A4C

### WP Figure 11C Ice industry

did_model <- feols(Pro_army ~ Business_Ice*treatment_factor | name  + treatment_factor,cluster = ~ name ,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)

coef_names

# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Business_Ice:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference  labs(title = "Difference-in-Differences Estimates Over Time",
  labs(title = "Ice maker",
       x = "",
       y = "") +
  custom_theme


print(plot)



```



```{r}

# Figure A4D

### WP Figure 11D Forestrty industry



did_model <- feols(Pro_army ~ Forestry*treatment_factor | name  + treatment_factor,cluster = ~ name ,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)

coef_names

# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Forestry:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference  labs(title = "Difference-in-Differences Estimates Over Time",
  labs(title = "Forestry",
       x = "",
       y = "") +
  custom_theme


print(plot)


```

```{r}

# Figure A5A

# WP Figure 12A: Bank (any)


did_model <- feols(Pro_army ~ Bank_Any*treatment_factor | treatment_factor + name,cluster = ~ name ,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)

# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Bank_Any:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference  labs(title = "Difference-in-Differences Estimates Over Time",
  labs(title = "Banks (Any)",
       x = "",
       y = "") +
  custom_theme

print(plot)



```

```{r}

# Figure A5B

### WP Figure 12B Mutual banks vs Other banks

BL <- c(as.Date('1937-03-22'),0,0,0,0)


did_model <- feols(Pro_army ~  Bank_Mutual*Bank_Any*treatment_factor  | treatment_factor+name,cluster = ~  name,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)


# Create a data frame for plotting
plot_data1 <- data.frame(
   Time = ifelse(grepl("^Bank_Mutual:treatment_factor", coef_names),
                          as.Date(gsub("Bank_Mutual:treatment_factor", "", coef_names), format = "%Y-%m-%d"),
                          as.Date(NA)),
#  Time = as.numeric(gsub("Textile_Wool_Linen:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)


# Calculate confidence intervals (95%)
plot_data1$CI_lower <- plot_data1$Estimate - 1.96 * plot_data1$SE
plot_data1$CI_upper <- plot_data1$Estimate + 1.96 * plot_data1$SE

plot_data1 <- rbind(plot_data1,BL)
plot_data1$Type <- 'Mutual Banks (Small)'

plot_data2 <- data.frame(
    Time = ifelse(grepl("^Bank_Any:treatment_factor", coef_names),
                          as.Date(gsub("Bank_Any:treatment_factor", "", coef_names), format = "%Y-%m-%d"),
                          as.Date(NA)),
  Estimate = coefficients,
  SE = se
)

plot_data2

# Calculate confidence intervals (95%)
plot_data2$CI_lower <- plot_data2$Estimate - 1.96 * plot_data2$SE
plot_data2$CI_upper <- plot_data2$Estimate + 1.96 * plot_data2$SE

plot_data2 <- rbind(plot_data2,BL)
plot_data2$Type <- 'Other Banks'



plot_data <- rbind(plot_data1,plot_data2)

library(ggplot2)

plot_data <- plot_data[!is.na(plot_data$Time), ]
plot_data$Time <- as.Date(plot_data$Time)

dodge_width <- 30

plot <- ggplot(plot_data, aes(x = Time, y = Estimate, color = Type, shape = Type)) +
  geom_point(position = position_dodge(width = dodge_width) )+                    # Points for the estimates
  geom_line(aes(linetype = Type),  position = position_dodge(width = dodge_width))+                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=0.8,   position = position_dodge(width = dodge_width)) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
  labs(title = "Mutual Banks (Small) vs Other Banks",
       x = "",
       y = "") +
  scale_color_manual(values = c("Mutual Banks (Small)" = "black", "Other Banks" = "darkgrey"), breaks = c("Mutual Banks (Small)", "Other Banks")) +   # Reverse order
#  scale_shape_manual(values = c("Procured" = 16, "Sanctioned" = 17), breaks = c("Sanctioned", "Procured")) +
#  scale_linetype_manual(values = c("Procured" = "solid", "Sanctioned" = "dashed"), breaks = c("Sanctioned", "Procured")) 
    custom_theme  +
  theme(legend.position = "bottom",
       # Enlarge legend title
    legend.text = element_text(size = 16),                  # Enlarge legend labels
    legend.key.size = unit(1.5, "lines"),                   # Increase legend key size
    legend.key.width = unit(1, "cm")                        # Adjust legend key width
  )
#scale_color_manual(values = c("blue", "green")) +  # Custom colors for cotton and linen
#scale_shape_manual(values = c(16, 17))             # Custom shapes for cotton and linen


print(plot)




```

```{r}

# Mutual Bank (comparison to normal bank is reported in Figure 18B)


did_model <- feols(Pro_army ~ Bank_Mutual*treatment_factor | treatment_factor + name,cluster = ~ name ,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)

# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Bank_Mutual:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference  labs(title = "Difference-in-Differences Estimates Over Time",
  labs(title = "Bank (mutual)",
       x = "",
       y = "") +
  custom_theme

print(plot)



```




```{r}

# Figure B1

### WP Figure 13 Procurement

BL <- c(as.Date('1937-03-22'),0,0,0,0)


did_model <- feols(Pro_army ~ Procured*treatment_factor | treatment_factor + name,cluster = ~ name,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)

# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Procured:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

#BL <- c(1936.03,0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-12-08'), linetype = "dashed", color = "black") +    # Zero line for reference
          geom_vline(xintercept = as.Date('1939-09-01'), linetype = "dashed", color = "black") +    # Zero line for reference
        geom_vline(xintercept = as.Date('1937-07-07'), linetype = "dashed", color = "gray") +    # Zero line for reference
  labs(title = "Firms in 1942 Procurement list: Estimates with Two-way FE",
       x = "Time Period",
       y = "Treatment Effect (Estimate)") +
  theme_minimal()


print(plot)




```

```{r}

# Figure B2

### WP Figure 24: Indsutry procured - Industry not procured


BL <- c(as.Date('1937-03-22'),0,0,0,0)


did_model <- feols(Pro_army ~  Procured_full*Industry_Any*treatment_factor  | treatment_factor+name,cluster = ~  name,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)


# Create a data frame for plotting
plot_data1 <- data.frame(
    Time = ifelse(grepl("^Industry_Any:treatment_factor", coef_names),
                          as.Date(gsub("Industry_Any:treatment_factor", "", coef_names),format = "%Y-%m-%d"),
                          as.Date(NA)),
  Estimate = coefficients,
  SE = se
)


# Calculate confidence intervals (95%)
plot_data1$CI_lower <- plot_data1$Estimate - 1.96 * plot_data1$SE
plot_data1$CI_upper <- plot_data1$Estimate + 1.96 * plot_data1$SE

plot_data1 <- rbind(plot_data1,BL)
plot_data1$Type <- 'Industry: Not in Procurement List'

plot_data2 <- data.frame(
        Time = ifelse(grepl("^Procured_full:Industry_Any:treatment_factor", coef_names),
                          as.Date(gsub("Procured_full:Industry_Any:treatment_factor", "", coef_names),format = "%Y-%m-%d"),
                          as.Date(NA)),
  Estimate = coefficients,
  SE = se
)


# Calculate confidence intervals (95%)
plot_data2$CI_lower <- plot_data2$Estimate - 1.96 * plot_data2$SE
plot_data2$CI_upper <- plot_data2$Estimate + 1.96 * plot_data2$SE
plot_data2 <- rbind(plot_data2,BL)
plot_data2$Type <- 'Industry: in Procurement List'



plot_data <- rbind(plot_data1,plot_data2)

library(ggplot2)

plot_data <- plot_data[!is.na(plot_data$Time), ]
plot_data$Time <- as.Date(plot_data$Time)

dodge_width=20

str(plot_data)



plot <- ggplot(plot_data, aes(x = Time, y = Estimate, color = Type, shape = Type)) +
  geom_point(position = position_dodge(width = dodge_width)) +                    # Points for the estimates
  geom_line(aes(linetype = Type),position = position_dodge(width = dodge_width)) +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2,position = position_dodge(width = dodge_width)) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +    # Zero line for reference +
  scale_color_grey() +
      geom_vline(xintercept = as.Date('1941-12-08'), linetype = "dashed", color = "darkgray") +    # Zero line for reference
          geom_vline(xintercept = as.Date('1939-09-01'), linetype = "dashed", color = "darkgray") +    # Zero line for reference
        geom_vline(xintercept = as.Date('1937-07-07'), linetype = "dashed", color = "gray") +    # Zero line for reference
  labs(title = "Industrial Firms: Procured in 1940-42 vs Not",
       x = "",
       y = "") +
  scale_x_date(date_labels = "%Y") + 
#  scale_color_manual(values = c("Procured" = "blue", "Sanctioned" = "green"), breaks = c("Sanctioned", "Procured")) +   # Reverse order
#  scale_shape_manual(values = c("Procured" = 16, "Sanctioned" = 17), breaks = c("Sanctioned", "Procured")) +
#  scale_linetype_manual(values = c("Procured" = "solid", "Sanctioned" = "dashed"), breaks = c("Sanctioned", "Procured")) 
        theme_minimal()  +
  theme(legend.position = "bottom",
       # Enlarge legend title
    legend.text = element_text(size = 16),                  # Enlarge legend labels
    legend.key.size = unit(1.5, "lines"),                   # Increase legend key size
    legend.key.width = unit(1, "cm")                        # Adjust legend key width
  )
#scale_color_manual(values = c("blue", "green")) +  # Custom colors for cotton and linen
#scale_shape_manual(values = c(16, 17))             # Custom shapes for cotton and linen

print(plot)




```





```{r}

# Figure B3A

### WP Figure 15A Automobile

did_model <- feols(Pro_army ~ Industry_Auto*treatment_factor | treatment_factor + name,cluster = ~  name,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)

# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Industry_Auto:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-12-08'), linetype = "dashed", color = "black") +    # Zero line for reference
          geom_vline(xintercept = as.Date('1939-09-01'), linetype = "dashed", color = "black") +    # Zero line for reference
        geom_vline(xintercept = as.Date('1937-07-07'), linetype = "dashed", color = "gray") +    # Zero line for reference
  labs(title = "Automobile",
       x = "",
       y = "") +
  custom_theme

print(plot)


```

```{r}

# Figure B3B

# WP Figure 15B Pharmacentical


did_model <- feols(Pro_army ~Pharma*treatment_factor | treatment_factor + name, cluster = ~ name, data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)

# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Pharma:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-12-08'), linetype = "dashed", color = "black") +    # Zero line for reference
          geom_vline(xintercept = as.Date('1939-09-01'), linetype = "dashed", color = "black") +    # Zero line for reference
        geom_vline(xintercept = as.Date('1937-07-07'), linetype = "dashed", color = "gray") +    # Zero line for reference
  labs(title = "Pharmaceutical",
       x = "",
       y = "") +
  custom_theme

print(plot)





```

```{r}

# Figure B3C

# WP Figure 15C Industrial Painting

did_model <- feols(Pro_army ~ Industry_Paint*treatment_factor | treatment_factor + name,cluster = ~  name,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)


# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Industry_Paint:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
    geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-12-08'), linetype = "dashed", color = "black") +    # Zero line for reference
          geom_vline(xintercept = as.Date('1939-09-01'), linetype = "dashed", color = "black") +    # Zero line for reference
        geom_vline(xintercept = as.Date('1937-07-07'), linetype = "dashed", color = "gray") +    # Zero line for reference
  labs(title = "Industrial Painting",
       x = "",
       y = "") +
  custom_theme


print(plot)




```

```{r}

# Figure B3D

### WP Figure 15D Ceramic and Glass


did_model <- feols(Pro_army ~ Business_Pottery*treatment_factor | treatment_factor + name,cluster = ~   name,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)


# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Business_Pottery:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-12-08'), linetype = "dashed", color = "black") +    # Zero line for reference
          geom_vline(xintercept = as.Date('1939-09-01'), linetype = "dashed", color = "black") +    # Zero line for reference
        geom_vline(xintercept = as.Date('1937-07-07'), linetype = "dashed", color = "gray") +    # Zero line for reference
#       geom_vline(xintercept = 1941.5, linetype = "dashed", color = "blue") +    # Zero line for reference
  labs(title = "Ceramic and Glass",
       x = "",
       y = "") +
   custom_theme

print(plot)





```



```{r}

# Figure B3E

### WP Figure 15E Shipbuilding

did_model <- feols(Pro_army ~ Industry_Ship_Misc*treatment_factor | treatment_factor + name,cluster = ~  name,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)

# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Industry_Ship_Misc:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
    geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-12-08'), linetype = "dashed", color = "black") +    # Zero line for reference
          geom_vline(xintercept = as.Date('1939-09-01'), linetype = "dashed", color = "black") +    # Zero line for reference
        geom_vline(xintercept = as.Date('1937-07-07'), linetype = "dashed", color = "gray") +    # Zero line for reference
  labs(title = "Shipbuilding",
       x = "",
       y = "") +
  custom_theme

print(plot)


```

```{r}

# Figure B3F

### WP Figure 15F Heavy Industry (Any)


did_model <- feols(Pro_army ~ Industry_Heavy*treatment_factor | treatment_factor + name,cluster = ~   name,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)


# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Industry_Heavy:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-12-08'), linetype = "dashed", color = "black") +    # Zero line for reference
          geom_vline(xintercept = as.Date('1939-09-01'), linetype = "dashed", color = "black") +    # Zero line for reference
        geom_vline(xintercept = as.Date('1937-07-07'), linetype = "dashed", color = "gray") +    # Zero line for reference
#       geom_vline(xintercept = 1941.5, linetype = "dashed", color = "blue") +    # Zero line for reference
  labs(title = "Heavy Industry (Any)",
       x = "",
       y = "") +
    custom_theme

print(plot)


```






```{r}

# Figure B4A

### WP Figure 16A: Livestock industry

table(ehd$Livestock)

did_model <- feols(Pro_army ~ Livestock*treatment_factor | name  + treatment_factor,cluster = ~ name ,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)

coef_names

# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Livestock:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-12-08'), linetype = "dashed", color = "black") +    # Zero line for reference
          geom_vline(xintercept = as.Date('1939-09-01'), linetype = "dashed", color = "black") +    # Zero line for reference
        geom_vline(xintercept = as.Date('1937-07-07'), linetype = "dashed", color = "gray") +    # Zero line for reference "Difference-in-Differences Estimates Over Time",
  labs(title = "Livestock",
       x = "",
       y = "") +
  custom_theme


print(plot)


```


```{r}

# Figure B4B

### WP Figure 16B: Fishery industry

did_model <- feols(Pro_army ~ Fisher*treatment_factor | name  + treatment_factor,cluster = ~ name ,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)

coef_names

# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Fisher:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-12-08'), linetype = "dashed", color = "black") +    # Zero line for reference
          geom_vline(xintercept = as.Date('1939-09-01'), linetype = "dashed", color = "black") +    # Zero line for reference
        geom_vline(xintercept = as.Date('1937-07-07'), linetype = "dashed", color = "gray") +    # Zero line for reference "Difference-in-Differences Estimates Over Time",
  labs(title = "Fishery",
       x = "",
       y = "") +
  custom_theme

print(plot)



```


```{r}


# Dairy industry: not reported (No effect, unstable pre-trend)

did_model <- feols(Pro_army ~ Food_Dairy*treatment_factor | treatment_factor + name,cluster = ~ name ,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)

# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Food_Dairy:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference  labs(title = "Difference-in-Differences Estimates Over Time",
  labs(title = "Dairy",
       x = "",
       y = "") +
  custom_theme

print(plot)



```




```{r}

# Figure 3

### WP Figure 17: Sanctioned vs Procured horse race

did_model <- feols(Pro_army ~  Procured_full*Sanctioned*treatment_factor | treatment_factor+name,cluster = ~  name,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)


# Create a data frame for plotting
plot_data1 <- data.frame(
    Time = ifelse(grepl("^Procured_full:treatment_factor", coef_names),
                          as.Date(gsub("Procured_full:treatment_factor", "", coef_names), format = "%Y-%m-%d"),
                          as.Date(NA)),
  Estimate = coefficients,
  SE = se
)


# Calculate confidence intervals (95%)
plot_data1$CI_lower <- plot_data1$Estimate - 1.96 * plot_data1$SE
plot_data1$CI_upper <- plot_data1$Estimate + 1.96 * plot_data1$SE
BL <- c(as.Date('1937-03-22'),0,0,0,0)
plot_data1 <- rbind(plot_data1,BL)
plot_data1$Type <- 'Procured'

plot_data2 <- data.frame(
      Time = ifelse(grepl("^Sanctioned:treatment_factor", coef_names),
                          as.Date(gsub("Sanctioned:treatment_factor", "", coef_names), format = "%Y-%m-%d"), as.Date(NA)),
  Estimate = coefficients,
  SE = se
)


# Calculate confidence intervals (95%)
plot_data2$CI_lower <- plot_data2$Estimate - 1.96 * plot_data2$SE
plot_data2$CI_upper <- plot_data2$Estimate + 1.96 * plot_data2$SE
BL <- c(as.Date('1937-03-22'),0,0,0,0)
plot_data2 <- rbind(plot_data2,BL)
plot_data2$Type <- 'Sanctioned'


plot_data3 <- data.frame(
        Time = ifelse(grepl("^Procured_full:Sanctioned:treatment_factor", coef_names),
                          as.Date(gsub("Procured_full:Sanctioned:treatment_factor", "", coef_names), format = "%Y-%m-%d"), as.Date(NA)),
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data3$CI_lower <- plot_data3$Estimate - 1.96 * plot_data3$SE
plot_data3$CI_upper <- plot_data3$Estimate + 1.96 * plot_data3$SE
BL <- c(as.Date('1937-03-22'),0,0,0,0)
plot_data3 <- rbind(plot_data3,BL)
plot_data3$Type <- 'Procured and Sanctioned'


plot_data <- rbind(plot_data1,plot_data2,plot_data3)

library(ggplot2)

plot_data <- plot_data[!is.na(plot_data$Time), ]
plot_data$Time <- as.Date(plot_data$Time)

dodge_width <- 30

plot <- ggplot(plot_data, aes(x = Time, y = Estimate, color = Type, shape = Type)) +
  geom_point(position = position_dodge(width = dodge_width)) +                    # Points for the estimates
  geom_line(aes(linetype = Type),position = position_dodge(width = dodge_width)) +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.3, 
                position = position_dodge(width = dodge_width)) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
  labs(title ="Sanctioned vs Procured vs Sanctioned&Procured",
       x = "",
       y = "") +
  scale_color_manual(values = c("Procured" = "darkgrey","Sanctioned" = "black","Procured and Sanctioned" = "gray")) +  # Light gray for both types
  scale_linetype_manual(values = c("Procured" = "dashed","Sanctioned" = "solid","Procured and Sanctioned" = "dotted")) + 
    custom_theme  +
  theme(legend.position = "bottom",
       # Enlarge legend title
    legend.text = element_text(size = 16),                  # Enlarge legend labels
    legend.key.size = unit(1.5, "lines"),                   # Increase legend key size
    legend.key.width = unit(1, "cm")                        # Adjust legend key width
  )

#scale_color_manual(values = c("blue", "green")) +  # Custom colors for cotton and linen
#scale_shape_manual(values = c(16, 17))             # Custom shapes for cotton and linen

print(plot)




```

```{r}

# Figure C1A

### WP Figure 18A Veteran's Association


did_model <- feols(Pro_army ~ Interest_Veteran*treatment_factor | treatment_factor + name,cluster = ~ name,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)


# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Interest_Veteran:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference  labs(title = "Difference-in-Differences Estimates Over Time",
  labs(title = "Veterans Association",
       x = "",
       y = "") +
  custom_theme

print(plot)




```




```{r}

# Figure C1B

# WPFigure 18B Colonial Interest Group


did_model <- feols(Pro_army ~ Interest_Colonial*treatment_factor | treatment_factor + name,cluster = ~ name,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)

# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Interest_Colonial:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference  labs(title = "Difference-in-Differences Estimates Over Time",
  labs(title = "Colonial Interest Group",
       x = "",
       y = "") +
  custom_theme



print(plot)


```



```{r}

# Figure C1C

### WP Figure 18C Legislators who studied abroad

did_model <- feols(Pro_army ~ study_abroad*treatment_factor | name  + treatment_factor,cluster = ~ name ,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)

coef_names

# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("study_abroad:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference  labs(title = "Difference-in-Differences Estimates Over Time",
  labs(title = "Studied Abroad",
       x = "",
       y = "") +
  custom_theme


print(plot)



```


```{r}

# Figure C1D

### WP Figure 18D University Graduates


did_model <- feols(Pro_army ~univ*treatment_factor | treatment_factor + name,cluster = ~ name,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)

# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("univ:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
  labs(title = "University Graduates",
       x = "",
       y = "") +
  custom_theme

print(plot)


```

```{r}

# Figure C1E

### WP Figure 18E  Bureaucrat

did_model <- feols(Pro_army ~ Ministry*treatment_factor | treatment_factor + name,cluster = ~ name,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)


# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Ministry:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference  labs(title = "Difference-in-Differences Estimates Over Time",
  labs(title = "Bureaucrat",
       x = "",
       y = "") +
  custom_theme

print(plot)




```

```{r}

# Figure C1F

### WP Figure 18F  Military Experience

did_model <- feols(Pro_army ~ Military*treatment_factor | treatment_factor + name,cluster = ~ name,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)


# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Military:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference  labs(title = "Difference-in-Differences Estimates Over Time",
  labs(title = "Military Experience",
       x = "",
       y = "") +
  custom_theme

print(plot)




```

```{r}

# Figure C1G

# WP Figure 18G Reporter


did_model <- feols(Pro_army ~ Reporter*treatment_factor | treatment_factor + name,cluster = ~ name,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)


# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Reporter:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference  labs(title = "Difference-in-Differences Estimates Over Time",
  labs(title = "Reporter",
       x = "",
       y = "") +
  custom_theme

print(plot)




```

```{r}

# Figure C1H

### WP Figure 18H Teacher / Lecturer


did_model <- feols(Pro_army ~ Educator*treatment_factor | treatment_factor + name,cluster = ~ name,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)


# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Educator:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference  labs(title = "Difference-in-Differences Estimates Over Time",
  labs(title = "Teacher / Lecturer",
       x = "",
       y = "") +
  custom_theme

print(plot)




```



```{r}

# Lawyer  Not reported in Paper (No d-i-d effect but pre trend is random)

table(ehd$Lawyer)

did_model <- feols(Pro_army ~ Lawyer*treatment_factor | name  + treatment_factor,cluster = ~ name ,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)

coef_names

# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Lawyer:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference  labs(title = "Difference-in-Differences Estimates Over Time",
  labs(title = "Lawyer",
       x = "",
       y = "") +
  custom_theme


print(plot)


```


```{r}

# Figure E1

### WP Figure 19 (Coupled with attached Geopackage file: Constituency Map made with QGIS)  : Constituency level data used to color the map
# Constituency level averages


library(dplyr)

# Calculate the average values of Sanctioned and Procured by dist_name (district)
average_values <- ehd %>%
  group_by(dist_name) %>%
  summarize(
    avg_sanctioned = mean(Sanctioned, na.rm = TRUE),
    avg_procured = mean(Procured_full, na.rm = TRUE),
    avg_Hayashi = mean(Hayashi, na.rm = TRUE),
    avg_Anti42 = mean(Anti42, na.rm = TRUE),
  )

# View the result
print(average_values)

write.csv(average_values, "par_district.csv", fileEncoding = "UTF-8")


```


```{r}

# Figure Q1A

### FIgure 20A Rubber


did_model <- feols(Pro_army ~ Rubber*treatment_factor | treatment_factor + name,cluster = ~   name,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)

coef_names

# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Rubber:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
#       geom_vline(xintercept = 1941.5, linetype = "dashed", color = "blue") +    # Zero line for reference
 labs(title = "Rubber",
       x = "",
       y = "") +
    custom_theme

print(plot)




```

```{r}

# Figure Q1B

### WP Figure 20B Cultivation Union


did_model <- feols(Pro_army ~ Cultivate*treatment_factor | treatment_factor + name,cluster = ~   name,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)

coef_names

# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Cultivate:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
#       geom_vline(xintercept = 1941.5, linetype = "dashed", color = "blue") +    # Zero line for reference
 labs(title = "Land Cultivattion Union (Crop Producers)",
       x = "",
       y = "") +
    custom_theme

print(plot)





```



```{r}

# Figure Q1C

### WP Figure 20C Agcicultural Association (Special)

did_model <- feols(Pro_army ~ Agri_Interest_Special*treatment_factor | treatment_factor + name, cluster = ~  name,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)

# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Agri_Interest_Special:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE


BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
  labs(title = "Agricultural Association: Specialized (non-rice)",
       x = "",
       y = "") +
  custom_theme


library(stargazer)
stargazer(did_model)

print(plot)


```

```{r}

# Figure Q1D

### WP Figure 20D Agcicultural Association (Generic): Likely Rice

did_model <- feols(Pro_army ~ Agri_interest_gen*treatment_factor | treatment_factor + name, cluster = ~  name,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)

# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Agri_interest_gen:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE


BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
  labs(title = "Agricultural Association: Generic (likely rice)",
       x = "",
       y = "") +
  custom_theme


library(stargazer)
stargazer(did_model)

print(plot)





```






```{r}

# Figure R1

### WP Figure 21: T-Test Results for IRAA Membership

t_test_results <- lapply(c("Sanctioned", "Procured",  "Infra_Any", 
                           "Industry_Any", "Transport_Any", "Service_Any", 
                           "Food_Any", "Newspaper_Any", "Finance_Any"), function(var) {
  t_test <- t.test(
    x = ehd$Yokusan41[ehd[[var]] == 1],
    y = ehd$Yokusan41[ehd[[var]] == 0],
    alternative = "two.sided"
  )
  cleaned_var <- gsub("_Any", "", var)
  data.frame(
    Variable = cleaned_var,
    t_statistic = t_test$statistic,
    p_value = t_test$p.value,
    mean_diff = diff(t_test$estimate)
  )
})


# Combine results into a single data frame
t_test_results_df <- do.call(rbind, t_test_results)

t_test_results_df$Variable
t_test_results_df$t_statistic



# Create a visualization using ggplot2
library(ggplot2)

plot <- ggplot(t_test_results_df, aes(x = reorder(Variable, t_statistic), y = t_statistic, fill = Variable)) +
  geom_bar(stat = "identity") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
  coord_flip() +
  labs(
#    title = "T-Test Results for Business x Pro-Army Factions: 1942 Election",
    x = "",
    y = "T-Statistic"
  ) +
    scale_fill_manual(values = c("lightgray","lightgray","lightgray","lightgray", "lightgray","darkgray","darkgray","lightgray","lightgray")) +
  theme_minimal() +
  theme(legend.position = "none")  # Removes the legend


print(plot)


```

```{r}

# Figure R2

### WP Figure 22: T-Test Results for Seat Retention in the 1942 Election


t_test_results <- lapply(c("Sanctioned", "Procured",  "Infra_Any", 
                           "Industry_Any", "Transport_Any", "Service_Any", 
                           "Food_Any", "Newspaper_Any", "Finance_Any"), function(var) {
  t_test <- t.test(
    x = ehd$Retain42[ehd[[var]] == 1],
    y = ehd$Retain42[ehd[[var]] == 0],
    alternative = "two.sided"
  )
  cleaned_var <- gsub("_Any", "", var)
  data.frame(
    Variable = cleaned_var,
    t_statistic = t_test$statistic,
    p_value = t_test$p.value,
    mean_diff = diff(t_test$estimate)
  )
})


# Combine results into a single data frame
t_test_results_df <- do.call(rbind, t_test_results)

t_test_results_df$Variable
t_test_results_df$t_statistic


# Create a visualization using ggplot2
library(ggplot2)

plot <- ggplot(t_test_results_df, aes(x = reorder(Variable, t_statistic), y = t_statistic, fill = Variable)) +
  geom_bar(stat = "identity") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
  coord_flip() +
  labs(
#    title = "T-Test Results for Business x Pro-Army Factions: 1942 Election",
    x = "",
    y = "T-Statistic"
  ) +
    scale_fill_manual(values = c("lightgray","lightgray","lightgray","lightgray", "lightgray","darkgray","darkgray","lightgray","lightgray")) +
  theme_minimal() +
  theme(legend.position = "none")  # Removes the legend


print(plot)


```

```{r}

# Figure V1

### WP Figure 23: T-Test: ”Proposal 85-6 on the Establishment of Imperial National Morality” 1944.2

# THose who were present on 1944.2 vote

eh9 <- ehd[!is.na(ehd$Endorse8506),]


t_test_results <- lapply(c("Sanctioned", "Procured", "Business_Any", "Infra_Any", 
                           "Industry_Any", "Transport_Any", "Service_Any", 
                           "Food_Any", "Newspaper_Any", "Finance_Any"), function(var) {
  t_test <- t.test(
    x = eh9$Endorse8506[eh9[[var]] == 1],
    y = eh9$Endorse8506[eh9[[var]] == 0],
    alternative = "two.sided"
  )
  cleaned_var <- gsub("_Any", "", var)
  data.frame(
    Variable = cleaned_var,
    t_statistic = t_test$statistic,
    p_value = t_test$p.value,
    mean_diff = diff(t_test$estimate)
  )
})


# Combine results into a single data frame
t_test_results_df <- do.call(rbind, t_test_results)

t_test_results_df$Variable
t_test_results_df$t_statistic



# Create a visualization using ggplot2
library(ggplot2)

plot <- ggplot(t_test_results_df, aes(x = reorder(Variable, t_statistic), y = t_statistic, fill = Variable)) +
  geom_bar(stat = "identity") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
  coord_flip() +
  labs(
#    title = "T-Test Results for Business x Pro-Army Factions: 1942 Election",
    x = "",
    y = "T-Statistic"
  ) +
    scale_fill_manual(values = c("lightgray","lightgray","lightgray","lightgray","lightgray", "lightgray","darkgray","darkgray","lightgray","lightgray")) +
  theme_minimal() +
  theme(legend.position = "none")  # Removes the legend


print(plot)



```

```{r}

# Figure V2

### WP Figure 24: T-Test: ”Proposal 85-9 on the Total Armament of the Population” 1944.2

eh9$Endorse8509 <- as.numeric(eh9$Endorse8509)

t_test_results <- lapply(c("Sanctioned", "Procured", "Business_Any", "Infra_Any", 
                           "Industry_Any", "Transport_Any", "Service_Any", 
                           "Food_Any", "Newspaper_Any", "Finance_Any"), function(var) {
  t_test <- t.test(
    x = eh9$Endorse8509[eh9[[var]] == 1],
    y = eh9$Endorse8509[eh9[[var]] == 0],
    alternative = "two.sided"
  )
  cleaned_var <- gsub("_Any", "", var)
  data.frame(
    Variable = cleaned_var,
    t_statistic = t_test$statistic,
    p_value = t_test$p.value,
    mean_diff = diff(t_test$estimate)
  )
})


# Combine results into a single data frame
t_test_results_df <- do.call(rbind, t_test_results)

t_test_results_df$Variable
t_test_results_df$t_statistic



# Create a visualization using ggplot2
library(ggplot2)

plot <- ggplot(t_test_results_df, aes(x = reorder(Variable, t_statistic), y = t_statistic, fill = Variable)) +
  geom_bar(stat = "identity") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
  coord_flip() +
  labs(
#    title = "T-Test Results for Business x Pro-Army Factions: 1942 Election",
    x = "",
    y = "T-Statistic"
  ) +
    scale_fill_manual(values = c("lightgray","lightgray","lightgray","lightgray","lightgray", "lightgray","darkgray","darkgray","lightgray","lightgray")) +
  theme_minimal() +
  theme(legend.position = "none")  # Removes the legend


print(plot)


```


```{r}


### Figure W1 / WP Figure 25 is in stock market section below

```




```{r}

# Figure W2A

### WP Figure 26A Liquor sector and nationalization shock (Placebo)

table(ehd$Liquor)

did_model <- feols(Pro_army ~ Liquor*treatment_factor | name  + treatment_factor,cluster = ~ name ,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)

coef_names

# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Liquor:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1937-3-31'), linetype = "dashed", color = "black") +    # Zero line for reference
 #     geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference  labs(title = "Difference-in-Differences Estimates Over Time",
  labs(title = "Liquor: Partially Nationalized April 1937",
       x = "",
       y = "") +
  custom_theme


print(plot)


```


```{r}

# Figure W2B

### WP Figure 26B Electricity industry and nationalization shock (placebo)

did_model <- feols(Pro_army ~ Infra_Electricity*treatment_factor | treatment_factor + name,cluster = ~ name,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)


# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Infra_Electricity:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1938-10-28'), linetype = "dashed", color = "black") +    # Zero line for reference
 #     geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference  labs(title = "Difference-in-Differences Estimates Over Time",
  labs(title = "Electricity(Partially Nationalized Oct 1938)",
       x = "",
       y = "") +
  custom_theme


print(plot)




```

```{r}

# Figure X1B

### WP Figure 27B Colonial Business

# (Figure X1A / WP Figure 27A is in stock section)


did_model <- feols(Pro_army ~ Business_Colonial*treatment_factor | treatment_factor + name,cluster = ~   name,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)


# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Business_Colonial:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
 #     geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
#       geom_vline(xintercept = 1941.5, linetype = "dashed", color = "blue") +    # Zero line for reference
  labs(title = "Colonial Businesses",
       x = "",
       y = "") +
   custom_theme

print(plot)





```

```{r}

# Figure X1C

### WP Figure 27C Manchurian Railway


did_model <- feols(Pro_army ~ Manchu_Railway*treatment_factor | treatment_factor + name,cluster = ~   name,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)


# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Manchu_Railway:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
 #     geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
#       geom_vline(xintercept = 1941.5, linetype = "dashed", color = "blue") +    # Zero line for reference
  labs(title = "Manchurian Railway",
       x = "",
       y = "") +
   custom_theme

print(plot)





```

```{r}

# Figure X1D

### WP Figure 27D: Manchurian Business (Excluding Manchurian Railway)


did_model <- feols(Pro_army ~ Manchu_Business*treatment_factor | treatment_factor + name,cluster = ~   name,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)


# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Manchu_Business:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

plot <- ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
 #     geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
#       geom_vline(xintercept = 1941.5, linetype = "dashed", color = "blue") +    # Zero line for reference
  labs(title = "Other Manchurian Business",
       x = "",
       y = "") +
   custom_theme

print(plot)





```




```{r}

# Businesses in Korea; not reported in paper (random pre trend, no effect, no prior)


did_model <- feols(Pro_army ~ Korea_Business*treatment_factor | treatment_factor + name,cluster = ~  name,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)


# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Korea_Business:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
  labs(title = "Business in Korea",
       x = "",
       y = "") +
  custom_theme

```



```{r}

# Sugar (Not presented in paper: small number of legislators (only 6); hit by sanction, but pre-trend is strong)


did_model <- feols(Pro_army ~ Food_Sugar_Misc*treatment_factor | name  + treatment_factor,cluster = ~ name ,  data = dfb)

# Summary of the model
summary(did_model)

summary_model <- summary(did_model)

coefficients <- summary_model$coefficients
se <- summary_model$se

coef_names <- names(coefficients)

coef_names

# Create a data frame for plotting
plot_data <- data.frame(
  Time = as.Date(gsub("Food_Sugar_Misc:treatment_factor", "", coef_names)),  # Extract period number from coefficient names
  Estimate = coefficients,
  SE = se
)

# Calculate confidence intervals (95%)
plot_data$CI_lower <- plot_data$Estimate - 1.96 * plot_data$SE
plot_data$CI_upper <- plot_data$Estimate + 1.96 * plot_data$SE

BL <- c(as.Date('1937-03-22'),0,0,0,0)

plot_data <- rbind(plot_data,BL)

library(ggplot2)

ggplot(plot_data, aes(x = Time, y = Estimate)) +
  geom_point() +                    # Points for the estimates
  geom_line() +                     # Line connecting the estimates
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2, size=1.4) + # Error bars for confidence intervals
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size=2.0) +    # Zero line for reference
    geom_vline(xintercept = as.Date('1940-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference
      geom_vline(xintercept = as.Date('1941-7-28'), linetype = "dashed", color = "black") +    # Zero line for reference  labs(title = "Difference-in-Differences Estimates Over Time",
  labs(title = "Sugar",
       x = "",
       y = "") +
  custom_theme

```



```{r}

# WP Table D1

### WP Table B1: Robustness Check 1

# Column1: Limiting dependent variables to Factions only (dff dataseet)

# Run the DiD regression
did_model_r1 <- feols(Pro_army ~ Treatment*Sanctioned  | treatment_period + name, data = dff2, cluster = ~name, fixef.rm = "infinite_coef")

# Summary of the model
summary(did_model_r1)

# Column 2: Limiting dependent variables to Parliamentary votes only (dfp dataset)

# Run the DiD regression
did_model_r2 <- feols(Pro_army ~ Treatment*Sanctioned  | treatment_period + name, data = dfp2, cluster = ~name, fixef.rm = "infinite_coef")

# Summary of the model
summary(did_model_r2)

# Column 3: Adding 2 extra events described in paper (dfg dataset)

did_model_r3 <- feols(Pro_army ~ Treatment*Sanctioned  | treatment_factor + name, data = dfg2, cluster = ~name, fixef.rm = "infinite_coef")

# Summary of the model
summary(did_model_r3)

# Column 4: Removing the events in 1940.9-1941.

did_model_r4 <- feols(Pro_army ~ Treatment*Sanctioned  | treatment_factor + name, data = dfg3B, cluster = ~name, fixef.rm = "infinite_coef")

# Summary of the model
summary(did_model_r4)


model_list <- list("Model 1" = did_model_r1, "Model 2" = did_model_r2, "Model 3" = did_model_r3, "Model 3" = did_model_r4)

# Output Table B1
modelsummary(model_list, output = "latex")

```

```{r}

# WP Table D2 

### WP Table B2 (Robustness check 2)


# Run the DiD regression (Non socialist only)
did_model_r5 <- feols(Pro_army ~ Treatment*Sanctioned  | treatment_factor + name, data = dfo2[dfo2$Socialist34==0,], cluster = ~name, fixef.rm = "infinite_coef")

# Summary of the model
summary(did_model_r5)

# Run the DiD regression (Non Tokyo legislators only)
did_model_r6 <- feols(Pro_army ~ Treatment*Sanctioned  | treatment_factor + name, data = dfo2[dfo2$prefecture!="東京",], cluster = ~name, fixef.rm = "infinite_coef")

# Summary of the model
summary(did_model_r6)

# Run the DiD regression (Non Tokyo legislators only)
did_model_r7 <- feols(Pro_army ~ Treatment*Sanctioned  | treatment_factor + name, data = dfo2[dfo2$Cabinet36==0,], cluster = ~name, fixef.rm = "infinite_coef")

# Summary of the model
summary(did_model_r7)

# Run the DiD regression
did_model_r8 <- feols(Pro_army ~ Treatment*Sanctioned  | treatment_factor + name, data = dfo2[dfo2$Business_Any==1,], cluster = ~name, fixef.rm = "infinite_coef")

# Summary of the model
summary(did_model_r8)

model_list <- list("Model 1" = did_model_r5, "Model 2" = did_model_r6, "Model 3" = did_model_r7, "Model 3" = did_model_r8)

# Output Table B2
modelsummary(model_list, output = "latex")


```




```{r}

table(eh1$Procured,  eh1$Seiyu34)
table(eh1$Sanctioned,  eh1$Seiyu34)

table(eh1$Procured,  eh1$Minsei34)
table(eh1$Sanctioned,  eh1$Minsei34)

```



```{r}

library(dplyr)

# Calculate percentage of 1s for each binary variable and get top 50
top_50 <- ehd %>%
  # Select numeric columns only
  dplyr::select(where(is.numeric)) %>%
  # Calculate percentage of 1s for each column
  summarise(across(everything(), ~ mean(. == 1) * 100)) %>%
  # Pivot longer to create variable name-value pairs
  pivot_longer(cols = everything(), names_to = "variable", values_to = "percentage_1s") %>%
  # Arrange in descending order of percentage
  arrange(desc(percentage_1s)) %>%
  # Select top 50
  slice_head(n = 500)

# Display results
print(top_50)

write.csv(top_50,"var3642.csv")




```


```{r}

# var3642 is manually cleaned 

```


```{r}

# Figure J1A

### WP Figure 1A

# Profile statistics(more than 20% of legislators)

major <- read.csv(here("major3645.csv"))

# Cleaned data of 

str(major)

plot <- ggplot(major, aes(x = reorder(Category, Percent), y = Percent)) +
  geom_bar(stat = "identity", fill = "lightgray") +
  coord_flip() + # Flip coordinates to make category names readable
  labs(title = "",
       x = "",
       y = "(%)") +
  theme_minimal() +
  theme(axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 16, hjust = 1))

print(plot)




```


```{r}

# Figure J1B

### WP Figure 1B

# Profile statistics(less than 20% of legislators)
career <- read.csv(here("career3645.csv"))

career1 <- career[career$Interest==1,]

plot <- ggplot(career1, aes(x = reorder(Category, Percent), y = Percent)) +
  geom_bar(stat = "identity", fill = "lightgray") +
  coord_flip() + # Flip coordinates to make category names readable
#  labs(title = "Background (>30 legislators:%)",
   labs(title = "",
       x = "",
       y = "(%)") +
  theme_minimal() +
  theme(axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 13, hjust = 1))

print(plot)



```


```{r}

# Figure J2

### WP Figure 2
# Business (sector) exposure of legislators

biz <- read.csv(here("biz3645.csv"))

dim(biz)

str(biz)

biz1 <- biz[biz$Percent>2.8,]
biz2 <- biz[biz$Percent<2.8,]

str(biz1)

plot1 <- ggplot(biz1, aes(x = reorder(Category, Percent), y = Percent)) +
  geom_bar(stat = "identity", fill = "lightgray") +
  coord_flip() + # Flip coordinates to make category names readable
#  labs(title = "Business Interest (>30 Legislators; %)",
  labs(title = "",
       x = "",
       y = "(%)") +
  theme_minimal() +
  theme(axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 10, hjust = 1))

plot2 <- ggplot(biz2, aes(x = reorder(Category, Percent), y = Percent)) +
  geom_bar(stat = "identity", fill = "lightgray") +
  coord_flip() + # Flip coordinates to make category names readable
# labs(title = "Business Interest (>10 Legislators; %)",
  labs(title = "",
       x = "",
       y = "(%)") +
  theme_minimal() +
  theme(axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 10, hjust = 1))

print(plot1)
print(plot2)


```






```{r}

# Figure F1

### WP Online Appendix C: Figure C1

# Balance test of roll call: 


library(tidyr)

# 77th to 79 sessions roll call

law79b <- read.csv(here("laws7779b.csv"))

# Limit the sample to those who were present in 1941 Feburuary session

eh4b <- ehd[ehd$B410223 ==1,]


# Making a data set of sanctioned and procured legislators

set <- as.data.frame(cbind(eh4b$Sanctioned , eh4b$Procured, eh4b$name))

colnames(set) <- c("Sanctioned","Procured","name")

data <- merge(law79b,set,by="name",all.x=TRUE)

data <- data[!is.na(data$Sanctioned) & !is.na(data$l1),]


# Perform t-tests and store results
results <- data.frame(
  Variable = character(),
  t_value = numeric(),
  p_value = numeric()
)


for (var in paste0("l", 1:10)) {
  t_test <- t.test(data[[var]] ~ data$Sanction)
  results <- rbind(results, data.frame(
    Variable = var,
    t_value = t_test$statistic,
    p_value = t_test$p.value
  ))
}

# Add significance level
results$Significance <- cut(
  results$p_value,
  breaks = c(-Inf, 0.001, 0.01, 0.05, 1),
  labels = c("***", "**", "*", "")
)

# Figure F1 A
# Visualize the results in a bar plot
library(ggplot2)

ggplot(results, aes(x = Variable, y = t_value, fill = Significance)) +
  geom_bar(stat = "identity", position = "dodge") +
  geom_text(aes(label = Significance), vjust = -0.5) +
  labs(
    title = "T-tests for Regislation Roll Call: Sanctioned",
    x = "",
    y = ""
  ) +
  theme_minimal() +
  scale_fill_grey(start = 0.8, end = 0.2) 



# Perform t-tests and store results
results <- data.frame(
  Variable = character(),
  t_value = numeric(),
  p_value = numeric()
)


for (var in paste0("l", 1:10)) {
  t_test <- t.test(data[[var]] ~ data$Procured)
  results <- rbind(results, data.frame(
    Variable = var,
    t_value = t_test$statistic,
    p_value = t_test$p.value
  ))
}

# Add significance level
results$Significance <- cut(
  results$p_value,
  breaks = c(-Inf, 0.001, 0.01, 0.05, 1),
  labels = c("***", "**", "*", "")
)

#Figure F1 B
# Visualize the results in a bar plot
library(ggplot2)

ggplot(results, aes(x = Variable, y = t_value, fill = Significance)) +
  geom_bar(stat = "identity", position = "dodge") +
  geom_text(aes(label = Significance), vjust = -0.5) +
  labs(
    title = "T-tests for Regislation Roll Call: Procured",
    x = "",
    y = ""
  ) +
  theme_minimal() +
  scale_fill_grey(start = 0.8, end = 0.2) 


# Figure C1 visualized below

```


```{r}

# Figure F2A

### Age balance: WP Figure C2A

ehd$birth <- as.numeric(ehd$birth)

ks_sanctioned <- ks.test(
  x = ehd$birth[ehd$Sanctioned == 1],  # "birth" for Sanctioned = 1
  y = ehd$birth[ehd$Sanctioned == 0]   # "birth" for Sanctioned = 0
)

# Procured KS Test
ks_procured <- ks.test(
  x = ehd$birth[ehd$Procured == 1],    # "birth" for Procured = 1
  y = ehd$birth[ehd$Procured == 0]     # "birth" for Procured = 0
)

# Print the results
print("KS Test for Sanctioned:")
print(ks_sanctioned)

print("KS Test for Procured:")
print(ks_procured)


library(ggplot2)

# Subset data
group1 <- ehd$birth[ehd$Sanctioned == 1]
group2 <- ehd$birth[ehd$Sanctioned == 0]

# Combine into a single data frame for ggplot
density_data <- data.frame(
  birth = c(group1, group2),
  group = factor(c(rep("Sanctioned = 1", length(group1)), 
                   rep("Sanctioned = 0", length(group2))))
)

# Plot density graph
plot <- ggplot(density_data, aes(x = birth, fill = group, color = group)) +
  geom_density(alpha = 0.3, adjust = 1.5) +  # Adjust smoothing with `adjust`
  labs(title = "Density Plot of Year of Birth by Sanctioned Status",
       x = "",
       y = "") +
  scale_fill_manual(values = c("Sanctioned = 1" = "darkgray", "Sanctioned = 0" = "lightgray")) +
  scale_color_manual(values = c("Sanctioned = 1" = "darkgray", "Sanctioned = 0" = "lightgray")) +
  custom_theme +
  theme(legend.position = "top")

print(plot)




```

```{r}

# Figure F2B

### Age balance: WP Figure C2B

# Subset data
group1 <- ehd$birth[ehd$Procured == 1]
group2 <- ehd$birth[ehd$Procured == 0]

# Combine into a single data frame for ggplot
density_data <- data.frame(
  birth = c(group1, group2),
  group = factor(c(rep("Procured = 1", length(group1)), 
                   rep("Procured = 0", length(group2))))
)

# Plot density graph
plot <- ggplot(density_data, aes(x = birth, fill = group, color = group)) +
  geom_density(alpha = 0.3, adjust = 1.5) +  # Adjust smoothing with `adjust`
  labs(title = "Density Plot of Year of Birth by Procurement Status",
       x = "",
       y = "") +
  scale_fill_manual(values = c("Procured = 1" = "gray", "Procured = 0" = "darkgray")) +
  scale_color_manual(values = c("Procured = 1" = "gray", "Procured = 0" = "darkgray")) +
  custom_theme +
  theme(legend.position = "top")

print(plot)

```

```{r}

# Online Appendix G

### WP Online Appendix D

# install.packages(c("tidyverse", "igraph", "reshape2"))

library(tidyverse)
library(igraph)
library(reshape2)

```


```{r}

# Load the dataset
data  <- read_csv("Q76.csv")
data2 <- read_csv("proposals76.csv")
data3 <- read_csv("Legislator_main.csv")

data3 <- dplyr::select(data3,c("name","Sanctioned","Procured","Seiyu34","Minsei34","Socialist34"))

data <- merge(data,data2,by="name")
data <- merge(data,data3,by="name")

summary(data)

# Extract the 19 dummy variables
dummy_vars <- data %>% dplyr::select(-c(name,Sanctioned,Procured,Seiyu34,Minsei34,Socialist34))

# Here we are checking if the two rows (names) share at least one '1' in any of the variables
adjacency_matrix <- (as.matrix(dummy_vars) %*% t(as.matrix(dummy_vars))) > 0

# Convert the adjacency matrix to a graph object
g <- graph_from_adjacency_matrix(adjacency_matrix, mode = "undirected", diag = FALSE)



```


```{r}

# Convert the adjacency matrix into a graph
g <- graph_from_adjacency_matrix(adjacency_matrix, mode = "undirected", diag = FALSE)

# Add names to the vertices
V(g)$name <- data$name
V(g)$party <- data$Seiyu34

# Define two colors for the binary party variable (0 = opposition, 1 = ruling)
party_colors <- c("black", "darkgray")  # Blue for opposition, Red for ruling

# Map the party variable (0, 1) to the corresponding colors
V(g)$color <- party_colors[V(g)$party + 1]  # Add +1 to match index (R uses 1-based indexing)

# summary()



```

```{r}

# Figure G1

# Graphs in WP Online Appendix D

### WP Figure D1: Former Seiyukai (Conservative) Legislators: Black

g <- delete.vertices(g, degree(g) == 0)

layout <- layout_with_fr(g)

# Plot the network
plot <- plot(g, 
     layout = layout,               # Use a different layout
     vertex.size = 2,                # Size of the vertices (dots)
     vertex.shape = "circle",        # Shape of the vertices
     vertex.label = NA,              # Hide the labels (names) from the plot for clarity
     edge.width = 1,                 # Thickness of the edges (lines)
     edge.color = "lightgray",            # Color of the edges
#     vertex.color = Sanctioned,          # Color of the vertices
     vertex.frame.color = NA,        # Label color for visibility
     vertex.color = adjustcolor(V(g)$color, alpha.f = 0.9),
 #     rescale = FALSE,               # Do not automatically rescale the layout
 #      xlim = range(layout[, 1]) * 1.2,   # Adjust x-axis limits to fill space
#     ylim = range(layout[, 2]) * 1.2,   # Adjust y-axis limits to fill space
     asp = 0)                       # Aspect ratio to stretch the plot

# Basic network metrics
cat("Number of vertices (names):", vcount(g), "\n")
cat("Number of edges (connections):", ecount(g), "\n")
cat("Is the graph connected?", is_connected(g), "\n")


print(plot)




```

```{r}

# Figure G2

### WP Figure D2: Former Minseito (Liberal) Legislators: Black

g <- graph_from_adjacency_matrix(adjacency_matrix, mode = "undirected", diag = FALSE)

# Add names to the vertices
V(g)$name <- data$name
V(g)$party <- data$Minsei34

# Define two colors for the binary party variable (0 = opposition, 1 = ruling)
party_colors <- c("black", "darkgray")  # Blue for opposition, Red for ruling

# Map the party variable (0, 1) to the corresponding colors
V(g)$color <- party_colors[V(g)$party + 1]  # Add +1 to match index (R uses 1-based indexing)

g <- delete.vertices(g, degree(g) == 0)

layout <- layout_with_fr(g)

# Plot the network
plot <- plot(g, 
     layout = layout,               # Use a different layout
     vertex.size = 2,                # Size of the vertices (dots)
     vertex.shape = "circle",        # Shape of the vertices
     vertex.label = NA,              # Hide the labels (names) from the plot for clarity
     edge.width = 1,                 # Thickness of the edges (lines)
     edge.color = "lightgray",            # Color of the edges
#     vertex.color = Sanctioned,          # Color of the vertices
     vertex.frame.color = NA,        # Label color for visibility
     vertex.color = adjustcolor(V(g)$color, alpha.f = 0.9),
 #     rescale = FALSE,               # Do not automatically rescale the layout
 #      xlim = range(layout[, 1]) * 1.2,   # Adjust x-axis limits to fill space
#     ylim = range(layout[, 2]) * 1.2,   # Adjust y-axis limits to fill space
     asp = 0)                       # Aspect ratio to stretch the plot

# Basic network metrics
cat("Number of vertices (names):", vcount(g), "\n")
cat("Number of edges (connections):", ecount(g), "\n")
cat("Is the graph connected?", is_connected(g), "\n")


print(plot)



```


```{r}

# Table H1

### WP Table E1: Lasso 1942

ehz <- read.csv("Legislator_main.csv")

ehz <- ehz %>%
  filter(!is.na(prefecture) & prefecture != "") %>%  # Remove empty rows
  mutate(id = row_number()) %>%                     # Add an ID to preserve rows
  pivot_wider(names_from = prefecture, 
              values_from = prefecture, 
              values_fn = length, 
              values_fill = 0)

eh <- ehz
eh <- eh[!(apply(eh, 1, function(x) all(is.na(x) | x == ""))), ]
eh <- eh[!(is.na(eh$Every42))&eh$EDRW42==1, ]

numeric_vars <- eh %>%
  select_if(is.numeric)

eh2 <- numeric_vars  %>%
  select_if(~ !any(is.na(.)))


include_vars <- c(
  
 "Otheragriculture_Misc", "Mashroom_Misc", "Agriinstitute_Other_Misc", "Agridevelop_Institute_Misc", "Agrifinance_Misc", "Tabacco_Misc",  "Chiken_Misc", "Fruits_Misc", "Charcoal_Misc", "Grain_Misc", "Tea_Misc",   "Agriland", "Cultivate", "Noukai", "Fisher", "Forestry", "Sericulture",   "Livestock", "Agri_Interest_Any", "Interest_Other_Misc",   "Interest_Wara_Misc", "Interest_Pharma_Misc", "Interest_Realestate_Misc",   "Interest_Finance", "Interest_Elec_Misc", "Interest_Salt_Misc",   "Interest_Food_Misc", "Interest_Plank",
  
 "Interest_Const", "Interest_Mining_Misc", "Interest_Tourism_Misc", "Interest_Infra_Misc",   "Interest_Publisher_Misc", "Interest_Employer", "Interest_Food", "Interest_Manuf",  "Interest_Commerce", "Interest_Transport",   "Cotton_Uni_Misc", "Silk_Uni_Misc", "Linen_Uni_Misc", "Wool_Uni_Misc",  "Layon_Uni_Misc", "Leather_Uni_Misc", "Silk_Any", "LiquorUni",   "Chukin_Misc", "Sankumi", "COOP", "CreditUni", "COC",  "Interest_Industry_Any", "Interest_Infra_Any",   "Road_League_Misc", "Train_Leagiue", "Const_League", 

"Union_Peasants", "Union_Labor", "Interest_Welfare", "Interest_Medical_Misc" ,
  
 "Oil_Sales", "Oil_Refine", "Oil_Any",   "Postoffice", "Ration", "Business_Ice", "Business_Pottery",   "Business_Colonial", "Procured_full",  "Sanctioned",  "Industry_Machine", "Industry_Auto",   "Industry_Metal", "Industry_Steel", "Industry_Heavy",   "Industry_Precision_Misc", "Industry_Electronics",   "Industry_Heavy_Machine_Misc", "Industry_Ship_Misc",  "Industry_Airplane_Misc", "Fertilizer",   "Industry_Chemical",  "Industry_Paint",  "Industry_Light_Pipe_Misc", "Industry_Light_Typewriter_Misc",   "Industry_Light_Can_Misc", "Industry_Light_Bicycle_Misc",   "Industry_Light_Watch_Misc",    "Rubber", "Industry_Any",   "Transport_Train", "Transport_Ship", "Transport_Land",   "Transport_Land_Warehouse", "Transport_Land_Truck",   "Transport_Bus", "Transport_Any",   "Mining_Energy", "Mining_Oil_Misc",  "Mining_Any",   "Pharma", 

  "Textile_Cotton", "Textile_Silk", "Textile_Silk_Likely",  "Textile_Weaver", "Textile_Spinner",   "Textile_Linen_Misc", "Textile_Wool_Misc", "Textile_Leather_Misc",   "Textile_Other_Misc",  "Textile_Any",   "Commerce_Wholesale", "Commerce_Wholesale_Trade",  "Commerce_Retail", "Commerce_Retail_Departmentstore",  "Commerce_Market", "Commerce_Merchant_Misc", "Commerce_Any", "Realestate", "Const_Cement", "Const_Const", "Const_Brick_Misc", "Const_Any",  "Service_Tourism", "Service_Hotel", "Service_Movie_Misc",   "Service_Onsen_Misc", "Service_Gamble", "Service_Any",   "Infra_Electricity", "Infra_Electricity_Hydro",   "Infra_Electricity_Supply", "Infra_Electricity_Lighting_Misc",   "Infra_Fuel_Gas", "Infra_Any",   "Bank_Any", "Bank_Small", "Bank_Trust_Misc",   "Bank_Commerce_Misc", "Bank_Agricultural", "Bank_Saving",   "Bank_Mutual", "Stock_Broker_Misc", "Stock_Exchange_Misc",   "Stock_Any", "Insurance_Any", "Insurance_Accident_Misc",   "Insurance_Life", "Finance_Any" ,
  
 "Food_Other_Rice_Misc", "Food_Other_Fruit_Misc",   "Food_Other_Wine_Misc", "Food_Other_Hokkaido_Misc", "Food_Silk",   "Food_Farmtool_Misc", "Food_Dairy", "Food_Tabacco_Misc",   "Food_Shoyu", "Food_Salt_Misc", "Food_Fish",   "Food_Starch_Misc", "Food_Flour_Misc", "Food_Sugar_Misc",   "Food_Any",  "Liquor",  "Wood_Paper", "Wood_Paper_Pulp_Misc", "Wood_Plank",  "Wood_Plank_Forest", "Wood_Plank_Sawing", "Wood_Plank_Charcoal_Misc",  "Wood_Any",
  
 "Newspaper_Any",   "Telecom_Any",    "Publisher_Any",     "Hospital", "Manchu_Business", "Manchu_Railway",   "Korea_Business", "Taiwan_Business_Misc",

 "samurai"	,	"adapt"	,	"inherit"	,	

  "todai", "kyodai", "keio", "waseda", "chuo", "meiji", "hosei",
  "nichidai", "kwansei", "cadet", "hitotsubashi", "Ritsumei_Misc", "Doshisha_Misc",
  "Toyo_Misc", "Senshu_Misc", "Univ_Tohoku_Misc", "Univ_Hokkaido_Misc", "tedai",
  "shidai", "national", "science", "medical", "literature", "phd", "grad",
  "us_degree", "us_study", "uk_study", "france_study_Misc", "germany_study",
  "asia_study_Misc", "study_other", "study_abroad",  "teaching_school",
  "agri_school", "commerce_school", "crafts_school", "soldier_school",
  "other_high", "lycee",
  
  "三重", "京都", "佐賀", "兵庫", "北海道", "千葉", "和歌山", "埼玉", "大分", "大阪",
  "奈良", "宮城", "宮崎", "富山", "山口", "山形", "山梨", "岐阜", "岡山", "岩手",
  "島根", "広島", "徳島", "愛媛", "愛知", "新潟", "東京", "栃木", "沖縄", "滋賀",
  "熊本", "石川", "神奈川", "福井", "福岡", "福島", "秋田", "群馬", "茨城", "長崎",
  "長野", "青森", "静岡", "香川", "高知", "鳥取", "鹿児島",
  
    "PC_三重県", "PC_京都府", "PC_佐賀県", "PC_兵庫県", "PC_北海道", "PC_千葉県",
  "PC_和歌山県", "PC_埼玉県", "PC_大分県", "PC_大阪府", "PC_奈良県", "PC_宮城県",
  "PC_宮崎県", "PC_富山県", "PC_山口県", "PC_山形県", "PC_山梨県", "PC_岐阜県",
  "PC_岡山県", "PC_岩手県", "PC_島根県", "PC_広島県", "PC_徳島県", "PC_愛媛県",
  "PC_愛知県", "PC_新潟県", "PC_東京府", "PC_栃木県", "PC_沖縄県", "PC_滋賀県",
  "PC_熊本県", "PC_石川県", "PC_神奈川県", "PC_福井県", "PC_福岡県", "PC_福島県",
  "PC_秋田県", "PC_群馬県", "PC_茨城県", "PC_長崎県", "PC_長野県", "PC_青森県",
  "PC_静岡県", "PC_香川県", "PC_高知県", "PC_鳥取県", "PC_鹿児島県",
    
    "CC_下関市", "CC_久留米市", "CC_京都市", "CC_今治市", "CC_仙台市", "CC_佐世保市",
  "CC_倉敷市", "CC_八代市", "CC_八幡市", "CC_八幡浜市", "CC_函館市", "CC_別府市",
  "CC_半田市", "CC_名古屋市", "CC_呉市", "CC_和歌山市", "CC_四日市市", "CC_堺市",
  "CC_塩竈市", "CC_大分市", "CC_大垣市", "CC_大津市", "CC_大牟田市", "CC_大阪市",
  "CC_奈良市", "CC_姫路市", "CC_宇和島市", "CC_宇治山田市", "CC_宇部市",
  "CC_宇都宮市", "CC_宮古市", "CC_宮崎市", "CC_宮津市", "CC_富山市", "CC_小樽市",
  "CC_小田原市", "CC_尾道市", "CC_山口市", "CC_岡谷市", "CC_岩国市", "CC_岸和田市",
  "CC_川内市", "CC_川口市", "CC_川崎市", "CC_布施市", "CC_帯広市", "CC_平市",
  "CC_広島市", "CC_延岡市", "CC_新宮市", "CC_新潟市", "CC_日立市", "CC_旭川市",
  "CC_明石市", "CC_札幌市", "CC_東京市", "CC_松山市", "CC_松本市", "CC_松阪市",
  "CC_横浜市", "CC_横須賀市", "CC_水戸市", "CC_沼津市", "CC_津山市", "CC_津市",
  "CC_浜松市", "CC_浦和市", "CC_熊本市", "CC_神戸市", "CC_福井市", "CC_福岡市",
  "CC_福島市", "CC_米沢市", "CC_若松市", "CC_諫早市", "CC_豊橋市", "CC_足利市",
  "CC_那覇市", "CC_都城市", "CC_金沢市", "CC_鎌倉市", "CC_長岡市", "CC_長崎市",
  "CC_長野市", "CC_静岡市", "CC_首里市", "CC_高松市", "CC_高田市", "CC_高知市",
  "CC_鳥取市", "CC_鶴岡市", "CC_鹿児島市", "CC_鹿屋市",
  
    "Military",   "Navy",  "General",  "Colonel",   "Lieutenant", 
  "Conscript",   "Accounting_Military",   "Clerical_military", 
  "Military_Army_Medic_Misc",  "Military_Army_Journalist_Misc",   "Military_Army_Logistics_Misc", "Military_Army_Engeneering_Misc",  "Military_Navy_Torpedo_Misc",  "Military_Army_Cannon",  "Military_Airforce_Misc",
  
  "Localgov_Tokyo", "Localgov_Osaka", "Localgov_Kyoto_Misc", "Localgov_Hokkaido", 
  "Interior",  "Governor", "Interior_Central", 
  "Police", "Localgov", 
  "Secretary", "Treasury", "Taxcollect", "Monopoly_Misc", 
  "Diplomat",  "Trainministry", 
  "Teleministry", 
  "Eduministry", 
   "Statministry", 
  "Indusministry", 
  "Otherministry_Audit_Misc", "Otherministry_Army_Misc", "Otherministry_Navy_Misc", 
  "Otherministry_Reconst_Misc", "Otherministry_Legal_Misc", "Otherministry_Colonial_Misc", 
  "Otherministry_Intel_Misc", "Ministry", "Ministry_Elite", "Ministry_Counsel",
  
  "Primary_Principal", "Primary_Teacher", "Primary_Staff", "Primary", "Educator", 
  "Middleschool_Teacher", "Specialschool_Teacher", "Highschool_Teacher", 
  "Colonialschool_Teacher_Misc", "Girlsschool_Teacher", "Agrischool_Teacher_Misc", 
  "Commerceschool_Teacher", "Secondary", 
  "Professor", "Lecturer", "Tertiary", 
  "Tertiary_Waseda", "Tertiary_Hitotsubashi_Misc", "Tertiary_Tokyo_Misc", 
  "Tertiary_Tsukuba_Misc", "Tertiary_Keio_Misc", "Tertiary_Nichidai", 
  "Tertiary_Meiji_Misc", "Tertiary_Rikkyo_Misc", "Tertiary_Chuo_Misc", 
  "Tertiary_Takushoku_Misc", "Tertiary_Kwansei_Misc", "Tertiary_Kangaku_Misc", 
  "Tertiary_Housei", "Tertiary_Senshu_Misc", "Tertiary_Toyo_Misc", 
  "Tertiary_Daito_Misc", "Tertiary_Taisho_Misc", "Tertiary_Doshisha_Misc",
  
    "Reporter", "Chiefreporter", "Editor", "Other_Journalist", "Journalist", 
  "Yomiuri", "Tokyo_Asahi_Misc", "Osaka_Asahi_Misc", "Tokyo_Nichinichi_Misc", 
  "Osaka_Mainichi_Misc", "Asahi", "Mainichi", "Jiji", "Hochi", "Yorozu_Misc", 
  "Yamato_Misc", "Niroku_Misc", "Chugai_Misc", "Dentsu_Misc", "Chuo_Shinbun", 
  "Miyako_Shinbun_Misc", "Kokumin_Shinbun_Misc", "Localnews", 
  "Tokyo_News", "Tabloid",
 
  "Lawyer", "Prosecutor", "Judge", 
  "Doctor", 
  "Blue_Collar", "White_Collar", "Zaibatsu_Worker", "Engeneer", 
  "Noneconomic_Job", "School_Worker", "Kindergarten_Misc", 
  "Public_Sector_Job", "Shugiin_Misc", "Nichigin_Misc", "NHK_Misc", 
  "Public_Sector_Local", "Public_Sector_National", "Legal_Clark", 
  "Religious_Job", "Religious_Job_Christian_Misc", 
  "Religious_Job_Budhist_Misc", "Religious_Jon_Shintoist_Misc", 
  "Welfare_Institution", "Accountant_Misc", 
  "Writer", "Thinktank", 
  "Pharmacist_Misc", "Dentist_Misc", "Veterinarian_Misc", 
  "Manchu_Admin", "Manchu_Career", "Manchu_Railway", 
  "Korea_Admin", "Korea_Career_Misc", 
  "Taiwan_Admin", "Taiwan_Career_Misc"
  
)

X <- as.matrix(eh2 %>% dplyr::select(dplyr::one_of(include_vars)))

threshold <- 0.01  # 1%

# Calculate the proportion of `1`s for each variable (column)
proportion_ones <- colMeans(X)  # Works because matrix stores 0/1 values as numeric

# Identify columns that meet the threshold
valid_columns <- proportion_ones >= threshold

# Filter the matrix to keep only valid columns
X <- X[, valid_columns]

y <- eh2$Every42

z <- eh2$Hayashi

V <- as.data.frame(cbind(X,y))


set.seed(20240725)

# Perform lasso regression using glmnet
lasso_model <- glmnet(X, y,  alpha = 1)

# Perform cross-validation to find the best lambda
cv_model <- cv.glmnet(X, y,  alpha = 1)

# Get the best lambda value
best_lambda <- cv_model$lambda.min

# Print the best lambda value
print(best_lambda)

lasso_coefficients <- coef(lasso_model, s = best_lambda)

# Print the coefficients
print(lasso_coefficients)

```

```{r}
# Table H2

### WP Table E2 lasso in 1937

ehc <- ehz
ehc <- ehc[!(apply(ehc, 1, function(x) all(is.na(x) | x == ""))), ]
ehc <- ehc[!(is.na(ehc$Hayashi))&(ehc$Diet37==1)|ehc$Diet36==1, ]

numeric_vars <- ehc %>%
  select_if(is.numeric)

na_summary <-numeric_vars  %>%
  summarize(across(everything(), ~ sum(is.na(.))))

print(na_summary)

numeric_vars <- numeric_vars[!is.na(numeric_vars$Elected36),]

eh3 <- numeric_vars  %>%
  select_if(~ !any(is.na(.)))



Xc <- as.matrix(eh3 %>% dplyr::select(dplyr::one_of(include_vars)))

threshold <- 0.01  # 1%

# Calculate the proportion of `1`s for each variable (column)
proportion_ones <- colMeans(Xc)  # Works because matrix stores 0/1 values as numeric

# Identify columns that meet the threshold
valid_columns <- proportion_ones >= threshold

# Filter the matrix to keep only valid columns
Xc <- Xc[, valid_columns]

# Every42, Incumbent42	,Former42 ,

yc <- eh3$Every42

zc <- eh3$Hayashi

Vc <- as.data.frame(cbind(Xc,zc))

set.seed(20240725)

# Perform lasso regression using glmnet
lasso_model <- glmnet(Xc, zc,  alpha = 1)
#lasso_model <- glmnet(X, y, family = "multinomial", alpha = 1)

# Perform cross-validation to find the best lambda
cv_model <- cv.glmnet(Xc, zc,    alpha = 1)
#cv_model <- cv.glmnet(X, y, family = "multinomial", alpha = 1)

# Get the best lambda value
best_lambda <- cv_model$lambda.min

# Print the best lambda value
print(best_lambda)

lasso_coefficients <- coef(lasso_model, s = best_lambda)

# Print the coefficients
print(lasso_coefficients)

#summary(lasso_coefficients)


```



```{r}

# Stock market data

```


```{r}


df <- read.csv(here("Stock3742.csv"))

```

```{r}

# Figure O1

### WP Figure 7


df <- df %>%
  mutate(Date = as.Date(paste(Year, Month, "28", sep = "-"))) # %>%
#  filter(Date >= as.Date("1940-02-01"))  # Filter the data to show only after 1940-02-01


# Define the date for the vertical line
vline_date1 <- as.Date("1937-07-07")
vline_date2 <- as.Date("1939-08-31")
vline_date3 <- as.Date("1941-12-07")
vline_date4 <- as.Date("1940-09-01")
vline_date5 <- as.Date("1941-07-01")
vline_date6 <- as.Date("1940-05-21")

# Plot the graph with multiple variables 
plot <- ggplot(df, aes(x = Date, y = 総指数S37)) +
  geom_line() +  # Set colors for Price and others
  geom_vline(xintercept = as.numeric(vline_date1), color = "gray", linetype = "dashed") +  # Add vertical line
    geom_vline(xintercept = as.numeric(vline_date2), color = "gray", linetype = "dashed") +  # Add vertical line
  geom_vline(xintercept = as.numeric(vline_date3), color = "gray", linetype = "dashed") +  # Add vertical line
    geom_vline(xintercept = as.numeric(vline_date4), color = "darkgray", linetype = "solid") +  # Add vertical line
    geom_vline(xintercept = as.numeric(vline_date5), color = "darkgray", linetype = "solid") +  # Add vertical line
      geom_vline(xintercept = as.numeric(vline_date6), color = "gray", linetype = "dashed") +  # Add vertical line
#      geom_vline(xintercept = as.numeric(vline_date6), color = "gray", linetype = "doshed") +  # Add vertical line
#   scale_color_manual(values = c("Price" = "blue", "Other" = "lightgray")) +  # Customize colors
  labs(x = "", y = "", title = "Tokyo Stock Exchange Stock Market Index (end-of-month)") +
  theme_minimal() +
  theme(legend.position = "none") 

print(plot)



```


```{r}

# Figure O2A

### WP Figure F1A Petrochemical


df2 <- df %>%
  mutate(Date = as.Date(paste(Year, Month, "28", sep = "-"))) #%>%
#  filter(Date >= as.Date("1940-02-01"))  # Filter the data to show only after 1940-02-01


# Reshape the data to long format, excluding Year, Month, and Date
df_long <- df2 %>%
  pivot_longer(cols = -c(Year, Month, Date), names_to = "Variable", values_to = "Value")

# Define the date for the vertical line
#vline_date <- as.Date("1940-07-28")
vline_date1 <- as.Date("1940-09-01")
vline_date2 <- as.Date("1941-07-01")
vline_date3 <- as.Date("1940-05-21")
vline_date4 <- as.Date("1939-09-01")
vline_date5 <- as.Date("1941-12-07")
vline_date6 <- as.Date("1938-08-31")
vline_date7 <- as.Date("1937-04-15")
vline_date8 <- as.Date("1940-11-30")
vline_date9 <- as.Date("1937-07-07")

# Plot the graph with multiple variables 
plot <- ggplot(df_long, aes(x = Date, y = Value, group = Variable)) +
  geom_line(aes(color = ifelse(Variable == "化學工業", "Price", "Other"), 
                alpha = ifelse(Variable == "化學工業", 1, 1)), size = 1)  +  # Set colors for Price and others
  geom_vline(xintercept = as.numeric(vline_date1), color = "darkgray", linetype = "dashed") +  # Add vertical line
    geom_vline(xintercept = as.numeric(vline_date2), color = "darkgray", linetype = "dashed") +  # Add vertical line
     geom_vline(xintercept = as.numeric(vline_date3), color = "lightgray", linetype = "dashed") +  # Add vertical line
   scale_color_manual(values = c("Price" = "black", "Other" = "lightgray")) +  # Customize colors
  labs(x = "", y = "", title = "Petrochemical") +
  theme_minimal() +
  theme(legend.position = "none") +
  theme(axis.title = element_text(size = 24))

plot <- ggplot(df_long, aes(x = Date, y = Value, group = Variable)) +
  # Plot gray lines first
  geom_line(data = df_long[df_long$Variable != "化學工業", ], 
            aes(color = "Other"), size = 1.2, linetype = "solid") +
  geom_line(data = df_long[df_long$Variable == "化學工業", ], 
            aes(color = "Price"), size = 1.5) +
  # Add vertical lines
  geom_vline(xintercept = as.numeric(vline_date1), color = "gray30", linetype = "dashed", size = 0.6) +
  geom_vline(xintercept = as.numeric(vline_date2), color = "gray30", linetype = "dashed", size = 0.6) +
  geom_vline(xintercept = as.numeric(vline_date3), color = "gray50", linetype = "dashed", size = 0.6) +
  # Customize colors
  scale_color_manual(values = c("Price" = "black", "Other" = "#d3d3d3")) +
  # Labels and theme adjustments
  labs(x = "", y = "", title = "Petrochemical") +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5, margin = margin(b = 10)),
    legend.position = "none",
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 10)
  )


print(plot)



```


```{r}

# Figure O2B

### WP Figure F1B Mining and Iron

plot <- ggplot(df_long, aes(x = Date, y = Value, group = Variable)) +
  # Plot gray lines first
  geom_line(data = df_long[df_long$Variable != "鑛業", ], 
            aes(color = "Other"), size = 1.2, linetype = "solid") +
  # Overlay black line for "鑛業"
  geom_line(data = df_long[df_long$Variable == "鑛業", ], 
            aes(color = "Price"), size = 1.5) +
  # Add vertical lines
  geom_vline(xintercept = as.numeric(vline_date1), color = "gray30", linetype = "dashed", size = 0.6) +
  geom_vline(xintercept = as.numeric(vline_date2), color = "gray30", linetype = "dashed", size = 0.6) +
  geom_vline(xintercept = as.numeric(vline_date3), color = "gray50", linetype = "dashed", size = 0.6) +
  # Customize colors
  scale_color_manual(values = c("Price" = "black", "Other" = "#d3d3d3")) +
  # Labels and theme adjustments
  labs(x = "", y = "", title = "Mining and Iron") +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5, margin = margin(b = 10)),
    legend.position = "none",
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 10)
  )


print(plot)



```

```{r}

# Figure O2C

### WP Figure F1C: Textile Weaving

plot <- ggplot(df_long, aes(x = Date, y = Value, group = Variable)) +
  # Plot gray lines first
  geom_line(data = df_long[df_long$Variable != "紡績", ], 
            aes(color = "Other"), size = 1.2, linetype = "solid") +
  # Overlay black line for "鑛業"
  geom_line(data = df_long[df_long$Variable == "紡績", ], 
            aes(color = "Price"), size = 1.5) +
  # Add vertical lines
  geom_vline(xintercept = as.numeric(vline_date1), color = "gray30", linetype = "dashed", size = 0.6) +
  geom_vline(xintercept = as.numeric(vline_date2), color = "gray30", linetype = "dashed", size = 0.6) +
  geom_vline(xintercept = as.numeric(vline_date3), color = "gray50", linetype = "dashed", size = 0.6) +
  # Customize colors
  scale_color_manual(values = c("Price" = "black", "Other" = "#d3d3d3")) +
  # Labels and theme adjustments
  labs(x = "", y = "", title = "Textile Weaving") +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5, margin = margin(b = 10)),
    legend.position = "none",
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 10)
  )


print(plot)


```

```{r}

# Figure O2D

### WP Figure F1D Synthetic Textile

plot <- ggplot(df_long, aes(x = Date, y = Value, group = Variable)) +
  # Plot gray lines first
  geom_line(data = df_long[df_long$Variable != "人絹人繊", ], 
            aes(color = "Other"), size = 1.2, linetype = "solid") +
  geom_line(data = df_long[df_long$Variable == "人絹人繊", ], 
            aes(color = "Price"), size = 1.5) +
  # Add vertical lines
  geom_vline(xintercept = as.numeric(vline_date1), color = "gray30", linetype = "dashed", size = 0.6) +
  geom_vline(xintercept = as.numeric(vline_date2), color = "gray30", linetype = "dashed", size = 0.6) +
  geom_vline(xintercept = as.numeric(vline_date3), color = "gray50", linetype = "dashed", size = 0.6) +
  # Customize colors
  scale_color_manual(values = c("Price" = "black", "Other" = "#d3d3d3")) +
  # Labels and theme adjustments
  labs(x = "", y = "", title = "Synthetic Textile") +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5, margin = margin(b = 10)),
    legend.position = "none",
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 10)
  )


print(plot)




```



```{r}
# Figure O3A

### WP Figure F2A Bank and Credit Union

plot <- ggplot(df_long, aes(x = Date, y = Value, group = Variable)) +
  # Plot gray lines first
  geom_line(data = df_long[df_long$Variable != "銀行信託", ], 
            aes(color = "Other"), size = 1.2, linetype = "solid") +
  geom_line(data = df_long[df_long$Variable == "銀行信託", ], 
            aes(color = "Price"), size = 1.5) +
  # Add vertical lines
  geom_vline(xintercept = as.numeric(vline_date1), color = "gray30", linetype = "dashed", size = 0.6) +
  geom_vline(xintercept = as.numeric(vline_date2), color = "gray30", linetype = "dashed", size = 0.6) +
  geom_vline(xintercept = as.numeric(vline_date3), color = "gray50", linetype = "dashed", size = 0.6) +
  # Customize colors
  scale_color_manual(values = c("Price" = "black", "Other" = "#d3d3d3")) +
  # Labels and theme adjustments
  labs(x = "", y = "", title = "Banks and Credit Union") +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5, margin = margin(b = 10)),
    legend.position = "none",
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 10)
  )


print(plot)
```

```{r}

# Figure O3B

### WP Figure F2B Insurance

plot <- ggplot(df_long, aes(x = Date, y = Value, group = Variable)) +
  # Plot gray lines first
  geom_line(data = df_long[df_long$Variable != "保險", ], 
            aes(color = "Other"), size = 1.2, linetype = "solid") +
  geom_line(data = df_long[df_long$Variable == "保險", ], 
            aes(color = "Price"), size = 1.5) +
  # Add vertical lines
  geom_vline(xintercept = as.numeric(vline_date1), color = "gray30", linetype = "dashed", size = 0.6) +
  geom_vline(xintercept = as.numeric(vline_date2), color = "gray30", linetype = "dashed", size = 0.6) +
  geom_vline(xintercept = as.numeric(vline_date3), color = "gray50", linetype = "dashed", size = 0.6) +
  # Customize colors
  scale_color_manual(values = c("Price" = "black", "Other" = "#d3d3d3")) +
  # Labels and theme adjustments
  labs(x = "", y = "", title = "Insurance") +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5, margin = margin(b = 10)),
    legend.position = "none",
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 10)
  )


print(plot)
```




```{r}

# Figure O3C

### WP Figure F2C Stock brokers (including stock exchange)

plot <- ggplot(df_long, aes(x = Date, y = Value, group = Variable)) +
  # Plot gray lines first
  geom_line(data = df_long[df_long$Variable != "取引所", ], 
            aes(color = "Other"), size = 1.2, linetype = "solid") +
  geom_line(data = df_long[df_long$Variable == "取引所", ], 
            aes(color = "Price"), size = 1.5) +
  # Add vertical lines
  geom_vline(xintercept = as.numeric(vline_date1), color = "gray30", linetype = "dashed", size = 0.6) +
  geom_vline(xintercept = as.numeric(vline_date2), color = "gray30", linetype = "dashed", size = 0.6) +
  geom_vline(xintercept = as.numeric(vline_date3), color = "gray50", linetype = "dashed", size = 0.6) +
  # Customize colors
  scale_color_manual(values = c("Price" = "black", "Other" = "#d3d3d3")) +
  # Labels and theme adjustments
  labs(x = "", y = "", title = "Stock brokers") +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5, margin = margin(b = 10)),
    legend.position = "none",
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 10)
  )


print(plot)



```


```{r}

# Figure O3D

### WP Figure F2D: Miscellaneous (Including International Trade)

plot <- ggplot(df_long, aes(x = Date, y = Value, group = Variable)) +
  # Plot gray lines first
  geom_line(data = df_long[df_long$Variable != "雑", ], 
            aes(color = "Other"), size = 1.2, linetype = "solid") +
  geom_line(data = df_long[df_long$Variable == "雑", ], 
            aes(color = "Price"), size = 1.5) +
  # Add vertical lines
  geom_vline(xintercept = as.numeric(vline_date1), color = "gray30", linetype = "dashed", size = 0.6) +
  geom_vline(xintercept = as.numeric(vline_date2), color = "gray30", linetype = "dashed", size = 0.6) +
  geom_vline(xintercept = as.numeric(vline_date3), color = "gray50", linetype = "dashed", size = 0.6) +
  # Customize colors
  scale_color_manual(values = c("Price" = "black", "Other" = "#d3d3d3")) +
  # Labels and theme adjustments
  labs(x = "", y = "", title = "Miscellaneous (Including International Trade)") +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5, margin = margin(b = 10)),
    legend.position = "none",
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 10)
  )




```







```{r}

# Paper and Pulp: Not reported in paper

plot <- ggplot(df_long, aes(x = Date, y = Value, group = Variable)) +
  # Plot gray lines first
  geom_line(data = df_long[df_long$Variable != "製紙", ], 
            aes(color = "Other"), size = 1.2, linetype = "solid") +
  geom_line(data = df_long[df_long$Variable == "製紙", ], 
            aes(color = "Price"), size = 1.5) +
  # Add vertical lines
#  geom_vline(xintercept = as.numeric(vline_date4), color = "gray30", linetype = "dashed", size = 0.6) +
#  geom_vline(xintercept = as.numeric(vline_date5), color = "gray30", linetype = "dashed", size = 0.6) +
  geom_vline(xintercept = as.numeric(vline_date1), color = "gray30", linetype = "dashed", size = 0.6) +
    geom_vline(xintercept = as.numeric(vline_date2), color = "gray30", linetype = "dashed", size = 0.6) +
  # Customize colors
  scale_color_manual(values = c("Price" = "black", "Other" = "#d3d3d3")) +
  # Labels and theme adjustments
  labs(x = "", y = "", title = "Paper") +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5, margin = margin(b = 10)),
    legend.position = "none",
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 10)
  )


print(plot)


```








```{r}


# Railroad: Not reported in paper

plot <- ggplot(df_long, aes(x = Date, y = Value, group = Variable)) +
  # Plot gray lines first
  geom_line(data = df_long[df_long$Variable != "鐵道電軌", ], 
            aes(color = "Other"), size = 1.2, linetype = "solid") +
  geom_line(data = df_long[df_long$Variable == "鐵道電軌", ], 
            aes(color = "Price"), size = 1.5) +
  # Add vertical lines
  geom_vline(xintercept = as.numeric(vline_date1), color = "gray30", linetype = "dashed", size = 0.6) +
  geom_vline(xintercept = as.numeric(vline_date2), color = "gray30", linetype = "dashed", size = 0.6) +
  geom_vline(xintercept = as.numeric(vline_date3), color = "gray70", linetype = "dashed", size = 0.6) +
  # Customize colors
  scale_color_manual(values = c("Price" = "black", "Other" = "#d3d3d3")) +
  # Labels and theme adjustments
  labs(x = "", y = "", title = "Railroad") +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5, margin = margin(b = 10)),
    legend.position = "none",
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 10)
  )


print(plot)

```






```{r}

# Figure O4A

### WP Figure F3A: Shipbuilding

plot <- ggplot(df_long, aes(x = Date, y = Value, group = Variable)) +
  # Plot gray lines first
  geom_line(data = df_long[df_long$Variable != "造船", ], 
            aes(color = "Other"), size = 1.2, linetype = "solid") +
  geom_line(data = df_long[df_long$Variable == "造船", ], 
            aes(color = "Price"), size = 1.5) +
  # Add vertical lines
  geom_vline(xintercept = as.numeric(vline_date4), color = "gray30", linetype = "dashed", size = 0.6) +
  geom_vline(xintercept = as.numeric(vline_date5), color = "gray30", linetype = "dashed", size = 0.6) +
  geom_vline(xintercept = as.numeric(vline_date1), color = "gray70", linetype = "dashed", size = 0.6) +
    geom_vline(xintercept = as.numeric(vline_date2), color = "gray70", linetype = "dashed", size = 0.6) +
      geom_vline(xintercept = as.numeric(vline_date9), color = "gray30", linetype = "dashed", size = 0.6) +
  # Customize colors
  scale_color_manual(values = c("Price" = "black", "Other" = "#d3d3d3")) +
  # Labels and theme adjustments
  labs(x = "", y = "", title = "Shipbuilding and Heavy Industry") +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5, margin = margin(b = 10)),
    legend.position = "none",
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 10)
  )


print(plot)



```


```{r}

# Figure O4B

### WP Figure F3B Machinery (Included Automobile)

plot <- ggplot(df_long, aes(x = Date, y = Value, group = Variable)) +
  # Plot gray lines first
  geom_line(data = df_long[df_long$Variable != "機械工作", ], 
            aes(color = "Other"), size = 1.2, linetype = "solid") +
  geom_line(data = df_long[df_long$Variable == "機械工作", ], 
            aes(color = "Price"), size = 1.5) +
  # Add vertical lines
  geom_vline(xintercept = as.numeric(vline_date4), color = "gray30", linetype = "dashed", size = 0.6) +
  geom_vline(xintercept = as.numeric(vline_date5), color = "gray30", linetype = "dashed", size = 0.6) +
  geom_vline(xintercept = as.numeric(vline_date1), color = "gray70", linetype = "dashed", size = 0.6) +
    geom_vline(xintercept = as.numeric(vline_date2), color = "gray70", linetype = "dashed", size = 0.6) +
      geom_vline(xintercept = as.numeric(vline_date9), color = "gray30", linetype = "dashed", size = 0.6) +
  # Customize colors
  scale_color_manual(values = c("Price" = "black", "Other" = "#d3d3d3")) +
  # Labels and theme adjustments
  labs(x = "", y = "", title = "Machinery and Automobile") +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5, margin = margin(b = 10)),
    legend.position = "none",
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 10)
  )


print(plot)



```

```{r}

# Figure O4C

### WP Figure F3C Ceramic and Glass

plot <- ggplot(df_long, aes(x = Date, y = Value, group = Variable)) +
  # Plot gray lines first
  geom_line(data = df_long[df_long$Variable != "窯業", ], 
            aes(color = "Other"), size = 1.2, linetype = "solid") +
  geom_line(data = df_long[df_long$Variable == "窯業", ], 
            aes(color = "Price"), size = 1.5) +
  # Add vertical lines
  geom_vline(xintercept = as.numeric(vline_date4), color = "gray30", linetype = "dashed", size = 0.6) +
  geom_vline(xintercept = as.numeric(vline_date5), color = "gray30", linetype = "dashed", size = 0.6) +
    geom_vline(xintercept = as.numeric(vline_date9), color = "gray30", linetype = "dashed", size = 0.6) +
  geom_vline(xintercept = as.numeric(vline_date1), color = "gray70", linetype = "dashed", size = 0.6) +
    geom_vline(xintercept = as.numeric(vline_date2), color = "gray70", linetype = "dashed", size = 0.6) +
  # Customize colors
  scale_color_manual(values = c("Price" = "black", "Other" = "#d3d3d3")) +
  # Labels and theme adjustments
  labs(x = "", y = "", title = "Ceramic and Glass") +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5, margin = margin(b = 10)),
    legend.position = "none",
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 10)
  )


print(plot)


```

```{r}

# Figure O4D

### WP Figure F3D Food

plot <- ggplot(df_long, aes(x = Date, y = Value, group = Variable)) +
  # Plot gray lines first
  geom_line(data = df_long[df_long$Variable != "食料品工業", ], 
            aes(color = "Other"), size = 1.2, linetype = "solid") +
  geom_line(data = df_long[df_long$Variable == "食料品工業", ], 
            aes(color = "Price"), size = 1.5) +
  # Add vertical lines
  geom_vline(xintercept = as.numeric(vline_date4), color = "gray30", linetype = "dashed", size = 0.6) +
  geom_vline(xintercept = as.numeric(vline_date5), color = "gray30", linetype = "dashed", size = 0.6) +
  geom_vline(xintercept = as.numeric(vline_date1), color = "gray70", linetype = "dashed", size = 0.6) +
    geom_vline(xintercept = as.numeric(vline_date2), color = "gray70", linetype = "dashed", size = 0.6) +
  # Customize colors
  scale_color_manual(values = c("Price" = "black", "Other" = "#d3d3d3")) +
  # Labels and theme adjustments
  labs(x = "", y = "", title = "Food") +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5, margin = margin(b = 10)),
    legend.position = "none",
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 10)
  )


print(plot)

```
```{r}

# Figure W1A

### WP Figure 25A: Stock Market Performance and Partial Nationalization Attempt (dashed line): General index in bold

plot <- ggplot(df_long, aes(x = Date, y = Value, group = Variable)) +
  # Plot gray lines first
  geom_line(data = df_long[df_long$Variable != "総指数S37", ], 
            aes(color = "Other"), size = 1.2, linetype = "solid") +
  geom_line(data = df_long[df_long$Variable == "総指数S37", ], 
            aes(color = "Price"), size = 1.5) +
  # Add vertical lines
  geom_vline(xintercept = as.numeric(vline_date8), color = "black", linetype = "dashed", size = 0.6) +
  # Customize colors
  scale_color_manual(values = c("Price" = "black", "Other" = "#d3d3d3")) +
  # Labels and theme adjustments
  labs(x = "", y = "", title = "Stock Index vs 'New Economic System' announcement") +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5, margin = margin(b = 10)),
    legend.position = "none",
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 10)
  )


print(plot)


```


```{r}

# Figure W1B

### WP Figure 25B: Electricity sector and Nationalization 


df2 <- df %>%
  mutate(Date = as.Date(paste(Year, Month, "28", sep = "-"))) #%>%
#  filter(Date >= as.Date("1940-02-01"))  # Filter the data to show only after 1940-02-01


# Reshape the data to long format, excluding Year, Month, and Date
df_long <- df2 %>%
  pivot_longer(cols = -c(Year, Month, Date), names_to = "Variable", values_to = "Value")

# Define the date for the vertical line
#vline_date <- as.Date("1940-07-28")
vline_date1 <- as.Date("1940-09-01")
vline_date2 <- as.Date("1941-07-01")
vline_date3 <- as.Date("1940-05-21")
vline_date4 <- as.Date("1939-09-01")
vline_date5 <- as.Date("1941-12-07")
vline_date6 <- as.Date("1938-08-31")
vline_date7 <- as.Date("1937-04-15")
vline_date8 <- as.Date("1940-11-30")
vline_date9 <- as.Date("1937-07-07")

plot <- ggplot(df_long, aes(x = Date, y = Value, group = Variable)) +
  # Plot gray lines first
  geom_line(data = df_long[df_long$Variable != "電燈電力", ], 
            aes(color = "Other"), size = 1.2, linetype = "solid") +
  geom_line(data = df_long[df_long$Variable == "電燈電力", ], 
            aes(color = "Price"), size = 1.5) +
  # Add vertical lines
  geom_vline(xintercept = as.numeric(vline_date6), color = "black", linetype = "dashed", size = 0.6) +
  # Customize colors
  scale_color_manual(values = c("Price" = "black", "Other" = "#d3d3d3")) +
  # Labels and theme adjustments
  labs(x = "", y = "", title = "Electricity") +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5, margin = margin(b = 10)),
    legend.position = "none",
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 10)
  )


print(plot)



```


```{r}

# Figure X1

### WP Figure 27A: Mamchurian Business

plot <- ggplot(df_long, aes(x = Date, y = Value, group = Variable)) +
  # Plot gray lines first
  geom_line(data = df_long[df_long$Variable != "満州", ], 
            aes(color = "Other"), size = 1.2, linetype = "solid") +
  geom_line(data = df_long[df_long$Variable == "満州", ], 
            aes(color = "Price"), size = 1.5) +
  # Add vertical lines
  geom_vline(xintercept = as.numeric(vline_date3), color = "black", linetype = "dashed", size = 0.6) +
  # Customize colors
  scale_color_manual(values = c("Price" = "black", "Other" = "#d3d3d3")) +
  # Labels and theme adjustments
  labs(x = "", y = "", title = "Manchhurian Businesses") +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5, margin = margin(b = 10)),
    legend.position = "none",
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 10)
  )


print(plot)


```





